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The  Analytic  Predictor  Corrector  (APC)  and  Energy  Controller  (EC)  atmospheric 
guidance  concepts  have  been  adapted  to  control  an  interplanetary  vehicle  aerobraking  in 
the  Martian  atmosphere.  Modifications  are  made  to  the  APC  to  improve  its  robustness  to 
density  variations.  These  modifications  include  adaptation  of  a  new  exit  phase  algorithm, 
an  adaptive  transition  velocity  to  initiate  the  exit  phase,  refinement  of  the  reference  dy¬ 
namic  pressure  calculation  and  two  improved  density  estimation  techniques.  The  modi¬ 
fied  controller  with  the  hybrid  density  estimation  technique  is  called  the  Mars  Hybrid 
Predictor  Corrector  (MHPC),  while  the  modified  controller  with  a  polynomial  density  esti¬ 
mator  is  called  the  Mars  Predictor  Corrector  (MPC). 

A  Lyapunov  Steepest  Descent  Controller  (LSDC)  is  adapted  to  control  the  vehicle. 
The  LSDC  lacked  robustness,  so  a  Lyapunov  tracking  exit  phase  algorithm  is  developed  to 
guide  the  vehicle  along  a  reference  trajectory.  The  equilibrium  glide  entry  phase  is  em¬ 
ployed  for  the  first  part  of  the  trajectory.  This  algorithm,  when  using  the  hybrid  density  es¬ 
timation  technique  to  define  the  reference  path,  is  called  the  Lyapunov  Hybrid  Tracking 
Controller  (LHTC).  With  the  polynomial  density  estimator  used  to  define  the  reference 
trajectory,  the  algorithm  is  called  the  Lyapunov  Tracking  Controller  (LTC). 


These  four  new  controllers  are  tested  using  a  six  degree  of  freedom  computer  simu¬ 
lation  to  evaluate  their  robustness.  MARS-GRAM  is  used  to  develop  realistic  atmo- 
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spheres  for  the  study.  The  atmospheres  are  then  perturbed  using  square  wave  density 
pulses.  The  MHPC,  MPC,  LHTC  and  LTC  show  dramatic  improvements  in  robustness 
over  the  APC  and  EC.  The  MHPC,  MPC,  LHTC  and  LTC  all  complete  the  initial  phase  of 
testing  (using  square  wave  density  pulses)  with  no  failures.  The  second  phase  tests  the 
MHPC,  MPC,  LHTC  and  LTC  against  atmospheres  where  the  inbound  and  outbound  den¬ 
sity  functions  are  different.  Square  wave  density  pulses  are  again  used,  but  only  for  the 
outbound  leg  of  the  trajectory.  Additionally,  sine  waves,  in  both  altitude  and  range,  are 
used  to  perturb  the  density  function.  All  four  new  controllers  are  able  to  compensate  for 
the  outbound  leg  density  pulses  with  no  hard  failures,  but  the  algorithms  are  sensitive  to 
large  amplitude  density  pulses.  Additionally,  these  control  algorithms  are  sensitive  to 
large  amplitude  sine  waves,  particularly  sine  waves  in  range.  The  hybrid  density  estimator 
responds  poorly  to  sine  waves  in  range  with  wavelength  between  twenty  and  two  hundred 
nautical  miles.  The  polynomial  density  estimator  is  sensitive  to  wavelengths  between  five 
hundred  and  two  thousand  nautical  miles.  Overall,  the  polynomial  density  estimator  per¬ 
forms  better  than  the  hybrid  density  estimator.  The  Lyapunov  tracking  phase  performs 
better  than  the  predictor  correctors  and  the  LTC  is  the  most  robust  control  algorithm  exam¬ 
ined. 
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CHAPTER  I 

INTRODUCTION 

When  orbital  transfer  is  required  near  a  celestial  body  with  an  atmosphere  of  suffi¬ 
cient  altitude  and  density,  it  is  often  advantageous  to  utilize  aerodynamic  forces  to  aid  in 
the  transfer*'*  ^  Aerodynamic  drag  forces  are  used  to  reduce  the  kinetic  energy,  while 
aerodynamic  lift  forces  are  used  to  control  the  trajectory  during  the  maneuver.  The  result 
is  a  vehicle  weight  savings  equivalent  to  the  propellant  necessary  to  perform  the  maneuver. 
The  critical  factor  for  success  in  the  aerobraking  maneuver  is  the  performance  of  the  guid¬ 
ance  control  system.  The  National  Aeronautics  and  Space  Administration  plans  a  1992 
launch  of  the  Aeroassisted  Flight  Experiment  (AFE)  to  serve  as  a  proof-of-concept  and 
test  vehicle  for  aerobraking  orbital  maneuvers*^.  Meanwhile,  an  aeroassisted  orbital 
transfer  maneuver  is  planned  for  the  Mars  Rover/Sample  Return  (MRSR)  Mission  to  re¬ 
duce  the  orbit  energy  from  the  hyperbolic  Martian  approach  orbit  to  capture  into  a  low 
Mars  orbit  with  a  commensurate  AV  savings  of  over  8000  ft/s  when  compared  with  an  all 
propulsive  mission  . 

Fig.  1  shows  the  proposed  interplanetary  mission  concept.  The  vehicle  will  be 
launched  from  Earth  into  an  elliptical  heliocentric  orbit.  The  vehicle  will  travel  almost 
eight  months  in  this  interplanetary  orbit  making  mid-course  corrections  as  necessary  to  in¬ 
tercept  Mars.  When  the  vehicle  reaches  Mars  there  will  be  6  km/sec  difference  between 
the  velocity  of  the  vehicle  and  Mars  orbital  velocity.  Without  some  method  of  changing 
the  vehicle’s  velocity  it  would  swing  by  Mars  without  capturing  into  a  Martian  orbit.  The 
proposed  method  for  imparting  this  velocity  change  is  to  use  the  aerodynamic  forces  im¬ 
parted  by  the  Martian  atmosphere.  The  final  mid-course  correction  to  the  interplanetary 
orbit  will  allow  the  vehicle  to  enter  the  Martian  atmosphere  as  shown  in  Fig.  2.  The  typi- 
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Fig.  1  Proposed  Interplanetary  Mission  Concept 
(Adapted  from  Reference  14) 
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cal  sequence  of  events  for  an  aerobraking  maneuver  call  for  the  vehicle  to  plunge  into  the 
atmosphere  and  fly  deep  in  the  atmosphere  until  the  velocity  is  appropriately  reduced. 
Then  the  vehicle  executes  a  pullout  maneuver  exiting  the  atmosphere  in  a  low  Mars  orbit. 
Finally  a  series  of  propulsive  maneuvers  are  performed  to  transition  the  vehicle  into  the 
desired  final  orbit. 


Fig.  2  Aerobraking  Maneuver  Sequence  of  Events 
(Adapted  from  Reference  14) 
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In  the  past,  for  space  missions  reentering  the  Earth’s  atmosphere,  only  the  destina¬ 
tion  coordinates  have  been  specified.  In  targeting  to  the  correct  orbit  following  the  aero- 
braking  maneuver  the  guidance  system  must  accurately  control  the  final  position  of  the 
vehicle  as  well  as  the  final  velocity  vector.  The  atmospheric  lift  and  drag  forces  affecting 
the  vehicle  are  proportional  to  the  atmospheric  density,  but  atmospheric  density  is  highly 
variable  .  The  guidance  algorithm  must  be  robust  enough  to  control  to  the  final  state 
even  with  these  uncertainties.  The  focus  of  this  dissertation  is  to  study  the  relative  merits 
of  several  existing  and  novel  guidance  algorithms,  with  particular  emphasis  upon  the  ex¬ 
tent  to  which  the  algorithms  tolerate  our  ignorance  of  the  Martian  atmosphere. 

Although  robustness  with  respect  to  density  variations  was  a  prime  factor  in  develop¬ 
ing  and  choosing  the  guidance  scheme  for  the  AFE^**®'*  becomes  even  more  criti¬ 

cal  for  the  Mars  mission.  Scientists  have  worked  for  many  years  to  characterize  the 
Earth’s  atmosphere.  Accelerometer  data  gathered  during  space  shuttle  returns  have  al¬ 
lowed  us  to  characterize  not  only  the  average  density  values  but  also  the  expected  magni¬ 
tude  and  frequency  of  the  random  density  variations'^’*^.  In  designing  the  AFE  guidance 
system  the  Earth’s  upper  atmosphere  in  the  region  250000-400000  feet,  was  assumed  to 
have  density  variations  of  ±25%  from  standard  values  over  small  altitude  intervals’^.  The 
Martian  atmosphere  goes  through  global  atmospheric  expansions  and  contractions  equiva¬ 
lent  to  an  atmospheric  shift  of  10  km^^.  Additionally,  data  gathered  from  the  Viking  1  and 
Viking  11  landers  show  density  variations  of  20  to  30%  over  small  altitude  intervals  in  the 
aerobraking  region’®’’^.  Since  we  have  only  two  sets  of  density  measurements  from  the 
Martian  aerobraking  region,  we  must  expect  even  larger  density  variations  to  occur.  It  is 
conjectured  that  density  shears  of  60%  or  greater  may  be  encountered.  Development  of  a 
robust  controller  capable  of  acceptable  performance  given  the  large,  unpredictable  density 
variations  in  the  Martian  atmosphere  is,  therefore,  vital  to  the  success  of  the  MRSR  mis¬ 


sion. 
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Background 

The  basic  technology  required  to  perform  hypersonic  flight  in  an  atmosphere  was  de¬ 
veloped  in  the  1950s  to  support  the  development  of  intercontinental  ballistic  missiles^^. 
Technology  was  extended  to  allow  the  Mercury  and  Gemini  projects  to  dissipate  kinetic 
energy  by  entering  the  atmosphere  with  low  ballistic  coefficient  vehicles.  Major  advances 
in  entry  technology  were  made  during  the  Apollo  program,  especially  in  the  areas  of  navi¬ 
gation,  guidance  and  control  during  atmospheric  maneuvering.  With  the  Space  Transpor¬ 
tation  System  (STS)  came  a  reusable  capability  to  deploy  and  retrieve  satellites  from  Low 
Earth  Orbit  (LEO).  Deployment  of  the  Space  Station  will  provide  a  permanent  base  in 
LEO  capable  of  performing  maintenance  and  repair  of  satellites.  However,  with  a  large 
percentage  of  satellites  in  Geosynchronous  Earth  Orbit  (GEO)an  economical  system  of 
deploying  satellites  to  GEO  and  then  returning  them  to  LEO  is  required.  The  National 
Aeronautics  and  Space  Administration  (NASA)  has  developed  an  aerobraking  vehicle,  the 
AFE,  to  meet  the  return  requirement^^.  In  designing  the  Mars  Rover/Sample  Return  mis¬ 
sion  an  aerobraking  phase  similar  to  that  of  the  AFE  is  envisioned  to  dissipate  kinetic  en¬ 
ergy  from  the  hyperbolic  Martian  approach  orbit  leaving  the  satellite  in  a  Low  Mars 
Orbit 

The  AFE,  scheduled  for  launch  in  1992  will  serve  as  a  proof  of  concept  and  test  ve¬ 
hicle  for  aerobraking.  AFE  will  enter  the  atmosphere  with  the  same  velocity  as  a  vehicle 
returning  from  GEO.  The  vehicle  will  fly  trimmed  at  a  constant  angle  of  attack,  and  there¬ 
fore,  at  near  a  constant  lift  to  drag  (L/D)  ratio.  AFE  will  roll  about  the  velocity  vector  to 
modulate  the  in  plane  portion  of  the  lift  to  control  the  trajectory  while  drag  dissipates  ki¬ 
netic  energy.  The  vehicle  will  exit  the  atmosphere,  after  an  appropriate  energy  reduction, 
with  exit  velocity  and  flight  path  angU  such  that  the  final  orbit  will  rendezvous,  at  apogee. 
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with  the  desired  LEO.  The  mission  concept  for  the  Mars  aerobrake  is  expected  to  be  quite 
similar,  but  of  course,  for  the  MRSR  objectives. 

The  complexity  of  the  competing  inequality  and  equality  constraints  placed  on  an 
aerobraking  maneuver  make  definition  of  an  optimal,  robust  control  algorithm  extremely 
difficult.  The  simplest  controllers  are  open  loop  controllers  designed  to  optimize  the  tra¬ 
jectory  for  a  specific  atmosphere,  entry  conditions  and  vehicle  design.  Talay,  et  al,^  opti¬ 
mized  a  bank  angle  history  for  a  nominal  1962  atmosphere  using  a  trajectory  optimization 
code.  When  this  bank  angle  history  was  used  in  trajectory  simulations  with  off-nominal 
atmospheres  in  several  cases  the  vehicle  either  exited  the  atmosphere  early  or  failed  to  exit 
at  all.  Vinh^  first  formulated  an  optimal,  minimum  fuel,  control  problem  using  a  com¬ 
bined  propulsive  and  aerodynamic  transfer.  He  shows  that  an  optimal  combined  propul- 
sive/aerodynamic  orbit  transfer  will  require  only  32%  of  the  total  AV  required  for  an  all 
propulsive  maneuver  for  an  orbit  transfer  from  GEO  to  LEO.  Then  Vinh,  et  al,^^  produce 
an  explicit  guidance  scheme  for  the  aerobraking  phase  of  a  drag  modulated  aeroassisted 
transfer  between  elliptical  orbits.  They  find  the  optimal  strategy  consists  of  bang-bang 
control  but  then  point  out  that  the  strategy  is  difficult  to  realize  because  the  switching  time 
must  be  very  accurate,  “within  a  fraction  of  a  second  to  avoid  crashing.”  They  propose  an 
alternative  strategy  whereby  the  drag  is  controlled  between  minimum  and  maximum  val¬ 
ues  as  a  function  of  the  current  state.  Kechichian,  et  al,^*  also  acknowledge  that  for  a  drag 
modulated  vehicle  bang-bang  control  is  optimum  for  minimizing  the  total  AV  required  to 
achieve  the  desired  orbit,  but  in  an  effort  to  reduce  the  sensitivity  to  switch  point  timing  a 
new  C£)n,ax’^Dmin*^Dmax  controller  is  developed  to  add  an  additional  degree  of  control. 
Sensitivity  analysis  shows  that  this  control  scheme  has  essentially  zero  sensitivity  to  an  at¬ 
mospheric  density  profile  of  ±15%  of  nominal  but  an  entry  corridor  width  of  ±0.1° 
should  be  maintained  to  avoid  excessive  AV  requirements. 
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Much  work  has  been  performed  in  the  area  of  optimal  aeroassisted  plane  changes. 
Hull,  Cl  al,  derives  an  optimal  guidance  scheme  for  performing  an  aeroassisted  plane 
change  between  circular  orbits.  They  assume  a  parabolic  drag  polar  for  the  vehicle  and 
use  Loh’s  con  tant^®  to  include  gravitational  terms  and  apparent  lift  terms  in  the  analysis. 
They  find  bank  angle  and  angle  of  attack  time  histories  which  minimize  the  total  AV  re¬ 
quired  to  perform  the  maneuver  by  maximizing  the  exit  velocity  following  the  aero  phase. 
Plane  changes  of  10  to  40  degrees  are  demonstrated.  Later  the  problem  is  reformulated^* 
using  heading  as  the  independent  variable  and  assuming  that  Loh’s  term  may  be  either 
positive  or  negative.  They  show  that  only  one  solution  exists  and  that  it  may  be  found  by 
solving  a  fourth  order  polynomial.  Hull,  McClendon,  and  Speyer^^  then  reform  the  prob¬ 
lem  assuming  an  elliptic  drag  polar  and  obtain  similar  results.  They  show  that  near  me 
end  of  the  atmospheric  turn  Loh’s  term  is  not  constant  which  may  cause  extremely  high 
angles  of  attack.  Finally,  Hull  and  co-workers^^  assume  Loh’s  term  is  piecewise  constant 
during  the  turn  and  reformulate  the  problem.  Using  the  method  of  successive  approxima¬ 
tions  they  construct  a  control  law  which  results  in  a  final  velocity  within  1%  of  the  true  op¬ 
timal  final  velocity  for  a  40°  plane  change  and  results  in  a  very  reasonable  maximum 
angle  of  attack  of  30°.  Johannesen,  et  al,^  formulate  an  approximate  control  h.w  for  lift 
and  bank  angle  to  maximize  orbit  plane  change  using  an  aeroassist  maneuver.  The  control 
law  is  tested  for  a  wide  range  of  speed  ratios  V^/V^.  They  observe  that  the  maximum  turn 
angle  for  any  speed  ratio  is  proportional  to  the  maximum  lift-to-drag  ratio. 

Two  unique  methods  of  determining  atmospheric  guidance  control  laws  have  been 
developed.  Mease  and  McCreary^  propo.se  using  an  approximate  closed  form  solution  of 
the  equations  of  motion.  Their  solution  divides  the  trajectory  into  three  regions.  During 
the  beginning  and  end  of  the  trajectory  the  gravitational  terms  are  assumed  to  dominate, 
while  in  the  mid-portion  of  the  trajectory  aerodynamic  terms  are  assumed  to  dominate. 
The  solutions  for  each  of  these  regions  is  combined  using  the  method  of  matched  asymp- 
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totic  expansion.  Final  apocenter  values  within  9%  of  the  targeted  values  are  demonstrated 
for  a  wide  range  of  entry  flight  path  angles  for  a  simulated  Mars  aerocapture  mission.  The 
other  unique  method  developed  by  Lee  and  Grantham^  uses  Lyapunov  optimal  feedback 
control  to  minimize  the  AV  required  to  raise  perigee  following  an  aerobraking  maneuver. 
This  method  calculates  a  descent  function  and  then  seeks  to  move  the  system  in  a  pre¬ 
ferred  direction,  opposite  the  gradient  of  the  descent  function.  The  Lyapunov  feedback 
controller  is  compared  with  an  optimal  open  loop  controller  derived  using  calculus  of  vari¬ 
ations  for  the  nominal  1962  standard  atmosphere.  Superior,  robust  performance  for  the 
Lyapunov  controller  is  demonstrated  for  both  the  standard  atmosphere  and  a  shuttle-de¬ 
rived  atmosphere. 

Control  laws  developed  using  optimal  control  theory  offer  excellent  perfonnance  in 
numerical  simulations,  but  those  methods  which  require  extensive  computation  for  each 
control  update  have  been  at  a  distinct  disadvantage  due  to  limitations  of  onboard  comput¬ 
ing  capability.  For  this  reason  several  simplified  guidance  schemes  have  been  developed. 
Letts  and  Pelekanos®  developed  a  control  law  using  bank-angle  modulation  of  the  lift  vec¬ 
tor  to  establish  a  constant  axial  deceleration  level  until  the  required  exit  velocity  is 
reached,  when  full  lift  up  is  commanded.  They  show  that  AV  required  to  circularize  fol¬ 
lowing  the  aerobrake  maneuver  increases  approximately  35  m/s  for  each  percent  change 
from  the  nominal  value  for  an  atmosphere  that  is  multiplied  by  a  constant  density  bias. 
Gamble,  et  al,^^  develop  a  control  .scheme  similar  to  that  of  Letts  and  Pelakanos  except 
Gamble’s  method  commands  to  an  equilibrium  glide  rather  than  a  constant  axial  decelera¬ 
tion  until  the  desired  velocity  is  achieved.  After  the  desired  velocity  is  achieved,  full  lift 
up  is  again  commanded  for  atmospheric  exit.  Ganible  finds  that  a  50%  increase  in  density 
has  little  effect  on  the  total  AV  required  but  a  .50%  decrease  in  density  increased  AV  re¬ 
quired  by  about  35%  due  to  problems  in  establishing  the  equilibrium  glide.  Cerimele,  et 
al,^*^  u.se  Gamble’s  equilibrium  glide  pha.se  during  the  entry  portion  of  the  trajectory,  but 
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then  switch  to  a  reference  drag  profile  like  Letts  and  Pelekanos  for  the  exit  portion  of  the 
trajectory.  Density  shears  in  the  atmosphere  are  simulated  and  the  A  V  required  following 
the  aerobraking  maneuver  is  found  to  be  very  sensitive  to  density  ratios  exceeding  ±30% 
occurring  over  altitude  ranges  of  1,000  to  10,000  ft.  Cerimele  and  Gamble^  produce  an 
analytic  predictor  corrector  guidance  algorithm,  again  using  the  equilibrium  glide  entry 
phase  but  with  a  predictor  corrector  exit  phase  designed  to  target  apogee  more  accurately. 
The  predictor  corrector  algorithm  assumes  a  constant  altitude  rate  and  an  exponential  at¬ 
mosphere  to  predict  apogee.  The  predictor  algorithm  iterates  altitude  rate  until  a  value  is 
found  which  produces  the  desired  apogee.  The  vehicle  is  then  commanded  to  this  altitude 
rate.  An  interesting  feature  added  to  this  algorithm  is  a  low  pass  density  filter.  Density  is 
computed  on-board  based  on  accelerometer  data.  Calculated  density  is  then  compared 
against  predicted  density  values  and  future  predicted  values  are  adjusted  accordingly.  This 
guidance  algorithm  was  tested  numerically  using  combined  dispersions  of  ±0.2°  in  entry 
flight  path  angle,  ±20%  density  variations  and  ±33%  L/D.  The  final  apogee  value  was 
within  2  nm  of  the  target  value  in  all  cases.  Gamble,  et  al,^^  present  three  atmospheric 
guidance  concepts  for  aeroassist  orbit  transfer  vehicles.  The  first  method  presented  is  the 
Analytic  Predictor  Corrector,  already  discussed.  The  second  is  a  Numerical  Predictor 
Corrector  algorithm  which  numerically  integrates  a  trajectory  assuming  constant  bank  an¬ 
gle  magnitude  and  an  assumed  density  profile.  The  bank  angle  is  iterated  until  the  desired 
apogee  is  computed  and  the  vehicle  is  commanded  to  this  bank  angle.  The  final  control  al¬ 
gorithm  presented  is  the  Energy  Controller  which  guides  the  vehicle  to  a  desired  energy 
state  at  atmospheric  exit.  The  energy  gain,  defined  as  the  ratio  of  energy  rate  to  energy  er¬ 
ror,  is  controlled  so  that  energy  error  exponentially  goes  to  zero  at  atmospheric  exit.  The 
energy  gain  command  is  converted  to  an  altitude  rate  command  which  in  turn  is  converted 
to  a  bank  angle  command.  Their  results  show  that  all  three  algorithms  are  capable  of 
maintaining  the  final  apogee  within  10  nm  and  AV  within  50  ft/sec  of  the  nominal  values 
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for  test  cases  with  dispersions  of  ±4  nm  in  perigee,  ±50%  in  density,  ±50%  in  W/Ci>\ 
and  +50%  in  L/D.  The  analytic  predictor  corrector  and  Energy  Controller  show  slightly 
worse  results  for  the  -50%  L/D  case. 

Fitzgerald  and  Ward*®''*  investigate  the  sensitivity  to  density  shears  of  the  Analytic 
Predictor  Corrector  and  Energy  Controller  algorithms  while  guiding  the  AFE  vehicle. 
They  consider  spike  and  step  shaped  density  dispersions  of  ±10  and  ±20%  magnitude 
with  durations  of  5,000  and  10,000  feet,  starting  at  altitudes  between  260K  and  295K  ft. 
A  V  increases  up  to  60%  for  the  Energy  Controller  and  41  %  for  the  APC  are  demonstrated. 
Fitzgerald**  then  formulates  a  Hybrid  Predictor  Corrector  algorithm  which  uses  the  atmo¬ 
spheric  density  profile  determined  during  the  entry  phase  in  the  predictor  corrector  of  the 
exit  phase.  This  significantly  reduces  the  sensitivity  to  density  shears  for  atmospheres 
where  the  exit  atmosphere  matches  the  entry  profile. 

Meyerson  and  Cerimele*^  review  the  aeroassist  vehicle  requirements  for  the  Mars 
Rover/Sample  Return  Mission.  They  use  a  modified  analytic  predictor  corrector  algorithm 
referred  to  as  HYPAS  as  the  controller  for  vehicles  with  L/D  ranging  from  0.3  to  1.5.  Ad¬ 
ditionally,  entry  velocities  from  5.79  to  9.20  km/sec  were  investigated.  A  recommenda¬ 
tion  of  this  study  is,  “to  refine  the  HYPAS  guidance  algorithm  to  control  the  trajectory 
more  accurately  in  the  exit  phase.”  They  recommend  using  two  exponential  atmosphere 
models  in  the  guidance  predictor. 

“The  Mars  atmosphere  is  highly  variable  on  a  daily,  seasonal  and  annual  basis*®.” 
The  thin  atmosphere  and  solar  heating  produce  a  large  daily  temperature  range  which 
translates  to  a  large  daily  density  fluctuation*®'^*.  Fig.  3  shows  that  at  the  surface  the 
Martian  atmospheric  density  is  approximately  two  orders  of  magnitude  less  than  that  of 
the  Earth’s  atmosphere  and  at  aerobraking  altitudes  there  is  still  more  than  an  order  of 
magnitude  difference  between  the  density  of  Earth  and  that  of  Mars.  “On  an  annual  basis. 
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Fig.  3  Comparison  of  COSPAR  Northern  Hemisphere  Mean  Mars 
Density  Model  and  1962  U.S.  Standard  Earth  Atmosphere 
(Adapted  from  Reference  19) 


spheric  pressure  at  the  surface  changes  by  ±15%  due  to  condensation  and  sublimation  of 
the  C02*^,”  which  produces  a  global  expansion  and  contraction  of  the  atmosphere  of 
roughly  10km.  Fig.  4  presents  the  density  deviations  of  the  COSPAR  high  density  model 
and  the  COSPAR  low  density  model  relative  to  the  COSPAR  Northern  Hemisphere  Mean 
Model.  Global  dust  storms  absorb  radiation  high  in  the  atmosphere,  thereby  increasing 
the  upper  atmosphere  temperature  and  causing  a  large  scale  expansion  of  the  atmosphere. 
The  density  is  then  substantially  increased  at  orbital  and  entry  altitudes.  Additionally, 
density  of  the  Martian  atmosphere  varies  widely  on  a  daily  basis.  Fig.  5  shows  the  expect¬ 
ed  morning  and  afternoon  density  profiles  calculated  for  summer  at  the  Viking  1  lander  lo¬ 
cation  while  Fig.  6  shows  the  calculated  density  profiles  for  winter  at  the  Viking  1  lander 
location.  These  figures  show  that  at  aerobraking  altitudes  the  density  may  vary  by  as 
much  as  100  to  150%  on  a  daily  basis.  The  Viking  1  and  2  landers  measured  atmospheric 
properties  during  their  descent  and  recorded  peak  to  peak  density  variations  in  the  aero- 


Fig.  4  Density  Deviations  of  the  COSPAR  Low-Cool  and  COSPAR 
High-Warm  Density  Models  as  Compared  to  the  COSPAR  Northern 

Hemisphere  Mean  Model 
(Adapted  from  Reference  19) 

tions  in  the  aerobraking  region  of  30%  over  a  15  km  altitude  band  and  20%  over  a  10  km 
region  *.  Fig.  7  and  Fig.  8  present  the  density  variations  measured  by  Viking  1  and  Vi¬ 
king  2  respectively  during  their  descent  to  Mars.  The  Mars  Global  Reference  Atmosphere 
Model  (MARS-GRAM)^^  characterizes  Mars  atmospheric  properties,  density,  tempera¬ 
ture,  pressure  and  wind  speed  and  direction,  as  functions  of  date,  time,  latitude,  longitude, 
altitude,  solar  activity  and  dust  storm  activity. 

This  dissertation  will  characterize  the  sensitivity  of  selected  aerobraking  guidance 
algorithms  with  respect  to  density  variations  of  the  type  and  magnitude  expected  in  the 
Martian  atmosphere  to  determine  their  suitability  to  perform  the  MRSR  Mission.  A  guid¬ 
ance  algorithm  capable  of  acceptable  performance  in  spite  of  the  uncertainty  in  Martian 
atmospheric  density  or  methods  of  reducing  the  uncertainty  will  be  developed.  To  attain 
this  goal  the  following  research  objectives  are  proposed. 


Fig.  5  Morning  and  Afternoon  Density  Profiles  Calculated  for  the 
Viking  1  Lander  Location  During  the  Summer 
(Adapted  from  Reference  19) 


Fig.  6  Morning  and  Afternoon  Density  Profiles  Calculated  for  the 
Viking  1  Lander  Location  During  the  Winter 
(Adapted  from  Reference  19) 
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Fig.  7  Viking  1  Measured  Density  Profile 
(Adapted  from  Reference  19) 


Fig.  8  Viking  2  Measured  Density  Profile 
(Adapted  from  Reference  19) 
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Research  Objectives 

When  designing  a  control  law  for  an  aeroassist  maneuver  an  exponential  variation  of 

atmosphere  density  with  altitude  is  an  extremely  attractive  computational  simplification. 

However,  given  the  large  density  biases  and  density  shears  of  the  Martian  atmosphere**' 
1  • 

,  a  guidance  algorithm  optimized  for  the  MRSR  aerobrake  flying  in  the  assumed  expo¬ 
nential  atmosphere  may  demonstrate  poor  performance  and  potentially  catastrophic  fail¬ 
ures  if  realistically  off-nominal  conditions  are  encountered*^.  Especially,  when  errors  in 
navigation  accuracy  and/or  vehicle  L/D  are  also  considered,  the  results  of  using  any  fixed 
density  model  may  be  catastrophic  with  the  vehicle  either  failing  to  be  captured  into  a 
Martian  orbit  or  failing  to  exit  the  Martian  atmosphere.  As  a  result,  the  sensitivity  of 
MRSR  atmospheric  guidance  to  perturbations  in  density  as  well  as  to  navigation  errors 
and  L/D  errors  is  critical  to  the  success  of  the  mission. 

A  systematic  method  of  evaluating  the  effects  of  density  biases  and  density  shears  in 
combination  with  navigation  errors  and  L/D  errors  on  an  MRSR  atmospheric  guidance  al¬ 
gorithm  is  sought.  The  methods  established  will  be  used  to  evaluate  candidate  guidance 
algorithms,  including  algorithms  developed  in  this  dissertation.  Toward  these  objectives, 
the  first  task  is  to  determine  the  expected  extremes  in  Martian  density.  MARS-GRAM^^ 
will  be  utilized  for  this  task.  The  highest  and  lowest  density  atmospheres  expected  are  de¬ 
termined  as  a  function  of  season,  time  of  day,  solar  activity  and  dust  storm  activity.  Since 
there  have  only  been  two  space  probes  which  have  measured  density  through  the  Martian 
aerobraking  region,  the  nature  and  magnitude  of  expected  worst  case  density  shears  is  not 
known  precisely  and  must  be  estimated.  These  atmosphere  extremes  are  checked  against 
the  Mars  standard  atmospheres**  and  the  Viking  data**'*^.  There  is  a  proposal  for  the 
Mars  Aeronomy  Observer  (MAO)^**'^*  to  send  additional  probes  to  the  Martian  surface  to 
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gather  more  data  quantifying  these  shears;  but,  this  mission  is  still  years  away.  The  ex¬ 
pected  navigational  accuracy  and  probable  L/D  errors  are  also  be  determined. 

Secondly,  the  Analytic  Predictor  Corrector  algorithm,^  the  algorithm  chosen  for  the 
AFE  mission,  and  the  Energy  Controller^'*  are  fine  tuned  for  the  MRSR  mission.  Then  the 
sensitivity  of  these  algorithms  when  faced  with  these  density  biases,  density  shears,  navi¬ 
gation  and  17D  errors  are  determined.  The  six  degree-of-freedom  Program  to  Optimize 
Simulated  Trajectories  (POST)^^  is  used  in  this  analysis.  The  sensitivity  of  these  algo¬ 
rithms  is  visually  presented  by  plotting  three  dimensional  sensitivity  surfaces  with  the 
qualitative  objective  of  finding  the  worst  combinations  of  dispersions  and  defining  the  per¬ 
formance  bounds  of  these  two  controllers.  Methods  of  improving  the  performance  of 
these  algorithms,  especially  methods  of  using  information  derived  early  in  the  trajectory, 
to  improve  the  performance  in  the  latter  portions  of  the  trajectory  (similar  to  the  methods 
proposed  by  Fitzgerald  in  his  Hybrid  Predictor  Corrector  algorithm)  are  evaluated.  Two 
new  algorithms  called  the  Mars  Hybrid  Predictor  Corrector  (MHPC)  and  the  Mars  Predic¬ 
tor  Corrector  (MPC)  are  developed.  These  two  algorithms  differ  only  in  their  density  esti¬ 
mation  techniques.  This  task,  along  with  developing  the  new  algorithm  proposed  below, 
are  crucial  to  the  research  and  secondary  only  to  the  task  of  defining  absolute  robustness 
limits. 

The  third  order  of  business  is  to  explore  more  elegant  ways  to  optimize  the  control¬ 
ler,  especially  ways  of  improving  the  robustness.  A  potential  candidate  Lyapunov  Steep¬ 
est  Descent  Controller^  (LSDC)  similar  to  the  one  suggested  by  Lee  and  Grantham  is 
coded  and  its  performance  tested  against  the  same  perturbations  as  the  others.  Two  new 
algorithms  are  developed,  again  differing  only  in  their  density  estimation  techniques. 
They  are  called  the  Lyapunov  Tracking  Controller  (LTC)  and  the  Lyapunov  Hybrid  Track¬ 
ing  Controller  (LHTC). 
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Finally,  and  most  importantly,  the  robustness  limits  of  the  improved  MFC,  MHPC, 
LTC  and  LHTC  controllers  are  characterized.  The  guidance  performance  is  thoroughly 
tested  to  find  the  tolerable  limits  on  density  bias  and  density  shear  given  the  probable  er¬ 
rors  in  navigation  and  L/D.  POST  is  used  to  test  the  guidance  algorithms,  using  the  Viking 
atmosphere  profiles**'*^.  The  limits  on  atmosphere  dispersions,  considering  the  inherent 
navigation  and  probable  L/D  errors,  under  which  acceptable  controller  performance  will 
occur  is  clearly  defined  from  the  results  of  these  simulations. 

These  limits  are  checked  against  the  worst  case  perturbations  expected  for  the  mis- 
sion**'*^’^^.  Conclusions  are  drawn  regarding  the  performance  of  these  algorithms  when 
faced  with  the  expected  density  variations,  as  well  as  possible  variations  in  vehicle  lift  to 
drag  ratio  and  entry  flight  path  angle.  Recommendations  for  future  study  are  then  presented. 

Organization  of  Dissertation 

Improvements  made  to  the  Analytic  Predictor  Corrector  and  Hybrid  Predictor  Cor¬ 
rector  control  algorithms  are  presented  in  Chapter  II.  Derivation  of  the  LSDC  and  a  LTA 
are  presented  in  Chapter  III.  Chapter  IV  details  the  model  used  for  the  trajectory  simula¬ 
tions  and  the  atmospheric  models.  Vehicle  mass  and  aerodynamic  data  are  presented 
along  with  atmospheric  entry  conditions.  The  controller  test  program  is  outlined  and  the 
perturbations  in  atmospheric  density,  vehicle  lift  and  drag,  and  entry  conditions  which 
were  used  in  the  test  program  are  presented.  Finally  the  performance  of  the  various  con¬ 
trollers  is  presented.  In  Chapter  V  the  four  best  performing  controllers,  the  MPC,  MHPC, 
LTA  and  the  LHTA  algorithms,  are  tested  against  each  other  in  a  head  to  head  fashion. 
The  perturbation  limits  which  define  the  edge  of  the  envelope  where  acceptable  perfor¬ 
mance  is  attainable  are  determined.  Conclusions  and  recommendations  are  presented  in 
Chapter  VI. 
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CHAPTER  II 

IMPROVEMENTS  TO  THE  PREDICTOR  CORRECTOR 

ALGORITHM 

The  Analytic  Predictor  Corrector  controller  (derived  in  Appendix  B)  is  the  control 
algorithm  selected  for  the  AFE*^.  This  controller  was  adapted  to  the  MRSR  problem  and 
tested  against  expected  perturbations  in  the  Martian  atmosphere  as  well  as  perturbations  in 
entry  conditions  and  vehicle  lift  and  drag  characteristics.  While  testing  the  Analytic  Pre¬ 
dictor  Corrector  controller,  several  areas  were  found  where  improvements  could  be  made. 
The  constant  multiplier  used  to  determine  the  reference  dynamic  pressure  was  changed  in 
an  effort  to  gain  robustness  and  prevent  premature  exit  from  the  atmosphere.  Borrowing  a 
technique  first  employed  by  Fitzgerald  an  improved  atmospheric  model  used  by  the  pre¬ 
dictor  step  to  determine  velocity  loss  during  the  exit  phase  was  also  incorporated.  Then  a 
new  method  of  estimating  density  incorporating  a  polynomial  to  fit  the  normalized  density 
function  was  developed^^.  A  modified  exit  phase,  first  developed  at  the  Charles  Stark 
Draper  Laboratory was  incorporated  and  tested,  first  without  and  then  with  the  im¬ 
proved  atmospheric  models.  The  new  exit  phase  also  assumes  a  constant  altitude  rate  to 
compute  the  velocity  loss  due  to  aerodynamic  drag.  However,  instead  of  predicting  the 
exit  state  and  iterating  altitude  rate  to  target  the  desired  apocenter,  the  velocity  required  to 
attain  the  desired  apocenter  altitude  is  computed  based  on  the  current  state  and  an  estimate 
of  the  remaining  velocity  loss  due  to  aerodynamic  drag.  The  iteration  is  simplified  to  a 
single  step  altitude  rate  correction.  With  the  large  density  shifts  present  in  the  Martian  at¬ 
mosphere  and  the  uncertainties  in  vehicle  and  entry  conditions  the  velocity  at  which  the 
controller  transitions  from  the  equilibrium  glide  phase  to  the  exit  phase  (incorporated  as  a 
controller  constant  for  the  ARC  controller  operating  in  the  Earth  atmosphere)  was  changed 
to  an  adaptive  parameter. 
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The  predictor  corrector  algorithm  developed  here  which  uses  a  variation  of  Fitzger¬ 
ald’s  density  estimation  scheme  will  be  referred  to  as  the  Mars  MHPC.  The  algorithm 
which  incorporates  the  polynomial  density  estimator  is  referred  to  simply  as  the  MPC  al¬ 
gorithm.  The  modifications  presented  here  convert  the  predictor  corrector  algorithm  from 
a  good  controller  for  guiding  the  aerobraking  phase  of  a  space  vehicle  returning  from  Geo¬ 
synchronous  Earth  Orbit  (GEO)  into  a  robust  control  algorithm  capable  of  accurately 
guiding  an  interplanetary  probe  through  an  aerobraking  maneuver  in  the  Martian  atmo¬ 
sphere. 


Reference  Dynamic  Pressure  Calculation 


The  equilibrium  glide  phase  of  the  APC  controller  seeks  an  equilibrium  condition 
with  the  vehicle  following  a  reference  dynamic  pressure  path.  The  reference  dynamic 
pressure  is  calculated  as  a  multiple,  K-,  of  the  dynamic  pressure  required  to  maintain 
equilibrium  with  the  lift  vector  oriented  down.  Gamble,  et  al^^,  recommend  that  this  mul¬ 
tiplier  be  1 .33  for  the  AFE  which  aerobrakes  in  the  Earth’s  atmosphere. 


^ref  = 
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Ideally,  to  have  the  minimum  AV  required  after  the  aerobraking  maneuver  the  veloc¬ 
ity  of  the  vehicle  should  be  decreased  as  high  in  the  atmosphere  as  possible.  Some  studies 
have  considered,  as  a  minimum  AV  aeromaneuver,  the  case  of  a  vehicle  with  infinite  lift 
skimming  the  edge  of  the  atmosphere  until  the  velocity  has  decreased  by  the  appropriate 
amount  so  the  vehicle  can  be  released  into  a  Hohman  Transfer  orbit  from  the  circular  orbit 
at  the  edge  of  the  atmosphere  to  the  desired  orbit^.  However,  decreasing  the  velocity  high 
in  the  atmosphere  means  flying  in  a  region  of  lower  density  and  consequently  lower  dy¬ 
namic  pressure.  Flying  higher  requires  the  in-plane  component  of  the  lift  vector  to  be  ori- 
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ented  more  downward  to  maintain  equilibrium.  This  idealization  is  satisfactory  only  with 
a  smooth  exponential  density  profile;  density  shears  are  not  allowed.  If  the  vehicle  is  fly¬ 
ing  in  equilibrium  using  50%  of  the  lift  down  capability  (O  =  ±120°)  and  a  density  shear 
is  encountered  which  decreases  the  density  by  50%  suddenly,  full  lift  down  will  be  re¬ 
quired  to  maintain  equilibrium.  In  actuality,  because  of  time  lags  between  the  encounter 
of  a  density  shear  and  the  vehicle’s  response,  coupled  with  the  potential  necessity  of  re¬ 
ducing  a  positive  alti  jde  rate,  the  minimum  acceptable  reference  dynamic  pressure  to 
maintain  control,  if  the  density  suddenly  decreases  by  50%,  is  considerably  more  than 
twice  that  required  to  maintain  equilibrium  in  a  full  lift  down  configuration. 

One  potential  drawback  to  increasing  K-  is  that  the  trajectory  loads  are  increased 
over  a  portion  of  the  flight.  Heat  rates  and  vehicle  acceleration  loads  are  increased  for  the 
portion  of  the  portion  of  the  trajectory  after  the  minimum  altitude  point  until  the  transition 
to  the  exit  phase,  however,  for  the  range  of  K-  between  1.33  and  4.5  the  maximum  heat 
rates,  g  loading,  the  minimum  altitude,  and  even  the  maximum  dynamic  pressure  do  not 
change.  4.5  was  the  largest  value  which  would  not  adversely  affect  the  peak  trajectory 
loads.  Furthermore,  it  has  been  found  that  the  total  heat  integrated  heat  load  calculated  is 
lower  for  a  higher  value  of  K-  because  the  vehicle’s  deceleration  is  greater  and  less  time  is 
required  to  reduce  the  vehicle’s  velocity. 

Fig.  9  presents  the  altitude  time  histories  and  Fig.  10  presents  the  velocity  time  histo¬ 
ries  for  trajectories  flown  through  a  nominal  Martian  atmosphere  perturbed  with  a  square 
wave  pulse  of  20,000ft  duration  located  between  140.000ft  and  160.000ft  altitude.  This 
pulse  multiplies  the  nominal  density  function  by  0.5  in  this  altitude  region.  The  three 
curves  presented  in  each  figure  represent  K~  values  of  1.33,  2.2  and  4.5.  Notice  that  the 

H 

trajectories  where  K-  equals  1.33  and  4.5  perform  well.  The  final  apocenter  is  within 
three  nautical  miles  of  the  targeted  270  nm  for  both  ca.ses.  However,  the  trajectory  flown 
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Fig.  9  Altitude  Time  Histories  Varied  K-  Values 


Fig.  10  Velocity  Time  Histories  Varied  K-  Values 
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with  K-  equal  to  2.2  skips  out,  or  exits  the  atmosphere  early  while  the  velocity  is  still  too 
high.  Skip  outs  are  difficult  to  predict  because  they  are  caused  by  a  sudden  negative  densi¬ 
ty  shear  which  reduces  the  vehicle’s  available  lift  and  drag  force.  If  this  negative  density 
shear  occurs  when  the  vehicle  is  in  a  relatively  .safe  regime  where  the  sudden  decrea.se  in 
density  will  not  place  the  vehicle  in  a  critical  situation  a  skip  out  may  not  occur.  However, 
if  the  vehicle  has  positive  altitude  rate,  or  is  performing  a  bank  reversal  and  the  lift  vector 
is  oriented  up,  or  the  control  scheme  allows  the  vehicle  to  overshoot  the  reference  dynamic 
pressure  trajectory,  thereby  temporarily  flying  at  a  dynamic  pressure  lower  than  the  refer¬ 
ence  value,  or  worst  yet,  combinations  of  the.se  factors  are  present  when  a  sudden  negative 
density  shear  is  encountered  a  skip  out  is  very  likely  to  occur.  Increasing  K~  will  not  al¬ 
ways  prevent  a  skip  out,  but  increasing  K-  does  tend  to  reduce  the  probability  of  a  skip 
out.  Indeed,  with  K-  set  to  4.5  (the  largest  value  possible  without  adversely  affecting  peak 
trajectory  loads)  no  skip  outs  were  encountered  during  the  test  program. 


With  the  uncertainties  in  the  Martian  atmosphere  the  slight  penalty  in  AV  required  to 
increase  this  multiplier  to  4.5  (less  than  10  ft/sec),  when  compared  to  the  1.33  value  rec¬ 
ommended  for  the  Earth  atmosphere,  seems  to  be  a  small  penalty  to  gain  additional  ro¬ 
bustness  and  limit  the  possibility  of  a  premature  .skip  out.  is  therefore  calculated 
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Improved  Exit  Phase  Density  Models 

The  second  area  of  improvement  is  the  density  estimation  technique.  Good  density 
estimation  is  critical  for  the  success  of  the  Martian  acrobrake.  With  the  wide  density 
swings  possible  in  the  Martian  atmosphere  the  correct  path  given  an  estimate  of  future 
density  may  prove  disastrous  if  that  estimate  is  wrong.  Given  that  the  MRSR  vehicle  will 


23 


traverse  over  1000  nautical  miles  during  the  aerobraking  maneuver  it  is  entirely  conceiv¬ 
able  that  the  density  function  encountered  may  vary  as  much  as  the  100  to  150%  variation 
in  density  between  morning  and  afternoon  presented  in  Fig.  5  and  Fig.  6.  The  density  es¬ 
timation  technique  must  not  only  develop  a  profile  of  density  versus  altitude,  but  must 
continue  to  update  this  estimate  throughout  the  trajectory  based  on  the  latest  accelerome¬ 
ter  based  density  measurement.  Two  methods  for  performing  this  task  are  presented  here. 

Hybrid  Density  Estimator 

A  method  similar  to  that  employed  by  Fitzgerald' '  to  model  the  atmosphere  is  em¬ 
ployed  here.  During  the  entry  phase  accelerometer-derived  density  is  recorded  near  each 
1000  ft  altitude  interval  along  with  the  altitude  for  each  density  measurement.  Then,  dur¬ 
ing  the  exit  phase  the  predictor  step  uses  these  measurements  to  predict  the  velocity  loss 
due  to  atmospheric  drag.  A  difference  between  this  approach  and  that  of  Fitzgerald  is  the 
inclusion  of  a  density  multiplier  derived  from  accelerometer-generated  density  measure¬ 
ments  which  continue  to  update  the  density  estimated  throughout  the  trajectory. 

The  accelerometer-generated  density  measurements  taken  during  the  entry  phase  of 
the  aerobraking  maneuver  is,  quite  likely,  the  best  estimate  of  the  atmospheric  density 
function  available  for  the  exit  phase  of  the  trajectory.  These  measurements  will  be  close, 
in  both  space  and  time,  to  the  density  for  the  remainder  of  the  flight.  They  will  indicate 
the  general  state  of  the  atmosphere,  that  is  whether  the  CO2  is  in  a  condensed  or  sublimat¬ 
ed  state,  and  they  will  provide  an  indication  of  the  vertical  wave  structure  of  the  atmo¬ 
sphere.  They  will  not  provide  any  information  on  horizontal  waves  which  may  affect 
density  on  the  outbound  leg.  To  compensate  for  this  latter  shortcoming  a  density  multipli¬ 
er  and  a  low  pass  filter  like  the  one  presented  for  the  APC^’  were  used.  However,  in- 
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stead  of  dividing  derived  density  by  the  density  predicted  by  an  exponential  function,  the 
divisor  is  the  density  predicted  from  the  measurements  taken  during  the  entry  phase. 

-(h-h.)/hS 

pRiodel  “ 

where  pj  is  the  density  which  was  measured  at  the  lower  edge  of  the  cunent  altitude  band, 
h]  is  the  altitude  at  which  this  measurement  was  taken,  and  hS  is  the  scale  height  for  the 
atmosphere  computed  between  this  density  measurement  and  the  next  measurement  ap¬ 


proximately  1000  ft  higher. 


■log(p2/pj)- 
.  *l“*2  . 


The  density  multiplier  is  then  computed  by  dividing  the  accelerometer  derived  densi¬ 
ty  by  the  density  predicted  for  the  current  altitude  using  the  density  model  derived  during 
entry.  The  result  is  filtered  using  a  low  pass  filter  to  remove  high  frequency  density  devia¬ 
tions  which  would  have  minimal  effect  on  the  post-aerobraking  apocenter. 


(l-/^)^o  +  ^(P/P„.odel) 


To  use  this  modified  atmosphere  in  the  predictor  step,  rewrite  Eq.  (148)  from  Appen¬ 
dix  B 

^  -(h-h^)/hSdh  ... 

—  =  -Cp.e  — .  (6) 


This  equation  may  be  integrated  assuming  a  constant  altitude  rate  to  give  the  velocity 
loss  due  to  atmospheric  drag  between  two  arbitrary  altitudes  h  j  and  /i2- 

r  -(hj-h,)/hS  -(A, -ti,)/hS, -]-l 

^r2  =  '  '  -e  '  '  (7) 

and,  with  the  scale  height  as  previously  calculated 
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r2 


Cpj  (h^-h^) 

/llog(p2/pj) 


(8) 


This  equation  gives  the  relative  velocity  at  /12  as  a  function  of  the  relative  velocity  at 
/i  j  and  the  densities  and  altitudes  at  the  two  locations.  The  method  for  employing  this  fea¬ 
ture  in  the  predictor  step  of  the  control  algorithm  is  to  first  use  the  velocity,  density  and  al¬ 
titude  at  the  current  satellite  location  as  the  subscript  1  variables  and  to  predict  the  velocity 
at  the  next  interval  where  density  measurements  were  stored  during  the  entry  using  that  al¬ 
titude  and  that  density  multiplied  by  the  density  multiplier  discussed  earlier  as  the  sub¬ 
script  2  variables.  Then  that  velocity  may  be  used  to  compute  the  velocity  at  the  next 
altitude  band  using  the  lower  stored  density  and  altitude  values  as  subscript  1  variables  and 
the  next  higher  density  and  altitude  measurements  as  subscript  two  variables.  Notice  that 
the  density  multiplier,  when  multiplied  by  each  of  the  stored  density  measurements,  will 
cancel  in  all  but  one  location. 


^^pP,  (^1  “^2^  ^P2  Y 

Alog(p2/pj)  Ui  J 


(9) 


This  procedure  is  repeated  until  the  exit  relative  velocity  is  computed.  The  velocity 
change  expected  between  the  current  location  and  atmospheric  exit  may  be  calculated  by 
subtracting  the  current  relative  velocity  from  the  predicted  exit  relative  velocity. 

(10) 


Polynomial  Density  Estimator 

The  second  method  of  density  estimation  curve  fits  a  sixth  order  polynomial  in  alti¬ 
tude  to  the  normalized  density  function.  This  technique  uses  accelerometer  derived  densi¬ 
ty  measurements  at  three  trajectory  locations  to  define  a  two  phase  exponential  function. 
Derived  density  is  recorded  at  one  second  intervals  and  then  normalized  by  the  exponen- 
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tial  function  in  an  effort  to  remove  the  underlying  predominant  exponential  component. 
As  the  satellite  reaches  the  bottom  of  the  trajectory  a  batch  estimate'^^  is  used  to  perform  a 
weighted  least  squares  fit  of  the  polynomial  coefficients  to  the  resulting  normalized  func¬ 
tion.  After  that,  a  sequential  estimate'^^  is  used  to  continue  updating  the  coefficients  of  the 
polynomial  for  the  remainder  of  the  trajectory. 

Based  on  MARS-GRAM  generated  data  a  two  phase  exponential  function  was 
chosen  to  normalize  the  density  data.  The  underlying  exponential  component  is  assumed 
to  be  two  exponential  functions  divided  at  250,000  ft  altitude  such  that  the  normalizing 
function  ^  is  expressed 

-  (A  -  250000)  /hS2|  ^  250, 000ft ) 

P(h)  =  *  '  .  (11) 

Pje  (A<250,000ft) 

hSl  and  hS2  are  the  scale  heights  below  and  above  250,000  ft.  is  the  density  at 
250,000  ft,  determined  using  accelerometer  derived  density  which  is  filtered  using  a  low 
pass  filter  like  the  one  presented  in  Eq.  (5).  The  scale  height  hS  1  is  found  by  using  the  fil¬ 
tered  density  measurement  when  the  vehicle’s  altitude  rate  first  becomes  positive  and  that 
at  250,000  ft  in  Eq.  (4).  similarly,  hS2  is  found  using  the  measured  density  at  400,000  ft 
and  the  measured  filtered  value  from  250,000  ft.  The  density  value  chosen  at  400,000  ft  is 
not  the  filtered  version  because  at  this  early  point  in  the  trajectory  the  density  filter  has  not 
had  sufficient  data  to  converge  to  a  reliable  estimate. 

After  the  altitude  rate  first  becomes  positive  and  the  constants  of  Eq.  (11)  have  been 
determined,  the  density  values  which  were  saved  at  one  second  intervals  during  the  de¬ 
scent  into  the  atmosphere  may  be  normalized.  The  resulting  data  is  fit  with  a  sixth  order 
polynomial  in  normalized  altitude  using  a  weighted  least  squares  (batch)  criterion  to  select 
the  coefficients  for  the  polynomial.  A  ninth  order  polynomial  was  originally  chosen  be- 
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cause  the  Viking  I  and  Viking  2  atmospheric  descent  data  (Fig.  7  and  Fig.  8)  shows  six 
and  five  major  density  extremes  respectively  in  the  aerobraking  region  and  the  Viking  1 
data  shows  four  additional  local  extremes.  This  would  require  at  least  a  seventh  order 
polynomial  to  model  even  the  major  extremes  accurately.  Because  computational  require¬ 
ments  for  the  density  filter  increase  approximately  as  the  square  of  the  order  of  the  polyno¬ 
mial  is  was  desired  to  use  as  low  order  as  practical.  After  some  testing  the  ninth  order 
polynomial  was  found  to  be  numerically  ill-conditioned.  A  sixth  order  polynomial  was 
found  to  be  much  better  behaved  and  adequate  for  modeling  the  expected  density  function. 
The  density  function  is  approximated  by 

p(/j)  =  p(h)  [Cj  C2J:  +  -I-  +  c^x^  +  Cf^x^  +  c^x^]  (12) 

h 

where  x  is  normalized  altitude,  jc  =  — . 

The  coefficients  Cj  through  c^  are  initially  determined  by  a  weighted  least  squares 
batch  estimate  presented  by  Junkins"*^.  The  procedure  is  to  begin  with  the  batch  of  m  nor¬ 
malized  density  measurements 


taken  at  the  m  known  altitude  locations  {hj)  at  m  one  second  intervals  until  the  altitude 
rate  becomes  positive.  The  altitudes  are  normalized  by  the  atmospheric  interface  altitude 
to  determine  the  x-s.  The  batch  estimator  must  select  the  coefficients  Cj  through  c-j  so 
that 


yy  =  c  1  -t-  C2XJ  +  c^xj  +  c^xj  +  c^Xj  +  c^Xj  +  c^Xj  +ej  (1 4) 

where  ej  is  the  residual  errors  after  selection  of  the  coefficients.  This  equation  may  be 
written  in  matrix  form 


Y  =  AC +  E 


(15) 
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where  t  is  the  estimate  of  the  coefficients  Cj  through  and  A  is  the  matrix 


-  1 

A 

■*1 

41 

A  = 

1 

^2 

A 

A 

4 

(16) 

1 

L 

4 

4 

An 

An 

4. 

The  batch  estimator  is  tasked  with  selecting  t  to  minimize  the  weighted  quadratic 
function  of  residual  errors 

(p  =  £^W'£-  (17) 

where  the  weighting  matrix  selected  is  a  diagonal  matrix  of  weights  applied  to  the  residual 
error  of  each  measurement. 


W  = 


Wj  0  0 

0  ^2  0 


L  0  0 


(18) 


Because  the  vehicle  is  traveling  into  the  region  of  higher  density  which  has  greater 
impact  on  the  sitell\e  trajectory  than  does  the  thinner  atmosphere  near  entry  and  exit  and 
because  more  recent  data  was  deemed  to  be  more  representative  of  future  density  than  was 
older  data,  the  weights  were  chosen  to  increase  with  time.  An  exponentially  increasing 
weighting  function  was  chosen  which  would  double  the  weight  after  1000  seconds.  This 
weighting  function  was  selected  through  experimentation  which  showed  that  a  slower  in¬ 
creasing  weighting  function  did  not  respond  quickly  enough  to  abrupt  density  shears  to 
produce  adequate  controller  performance,  while  faster  increasing  weighting  functions 
tended  to  ignore  data  gathered  early  in  the  trajectory  and  produced  a  poor  estimate  of  den¬ 
sity  in  the  upper  altitude  regions  which  also  had  a  negative  impact  on  controller  perfor¬ 
mance.  The  weights  selected  were 
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6.9314xl0~^/, 

w.  =  e  \  (19) 

Equation  ( 1 5)  is  solved  for  E  and  the  result  substituted  into  Eq.  ( 1 7).  After  some  ma¬ 
nipulation  ({>  is  expressed 

iSf  =  A^WAt.  (20) 

To  minimize  <j>  it  is  necessary  that 

V^(p  =  -2A'^H'y'-t-2A'’'vV'A^  =  0.  (21) 

This  equation  is  solved  to  obtain  the  weighted  least  squares  normal  equation  for  t 

t  =  {A^WA)~^A^WY.  (22) 

After  the  satellite  altitude  rate  becomes  positive  the  estimator  switches  to  a  linear  se¬ 
quential  estimator'*^.  To  facilitate  this  switch,  the  covariance  matrix  P  is  recorded  from 
the  batch  estimate 

P,  =  (23) 

where  for  the  first  sequential  estimation  step  the  k  subscripts  are  simply  the  matrix  values 
from  the  batch  estimate.  For  subsequent  steps  the  k  subscripts  will  indicate  values  from 
the  previous  step  while  k+l  will  indicate  updated  values.  A  linear  Kalman  filter  is  then 
employed  to  update  the  estimates  of  As  new  density  measurements  are  made  available 
at  one  second  intervals  the  estimate  of  t  is  updated  using 

For  the  sequential  estimator  W  is  just  the  scalar  value  of  w-  calculated  using  Eq.  (19).  A  is 
only  the  new  row  of  the  A  matrix  shown  in  Eq.  (16)  calculated  using  the  current  value  of  h. 
Y  is  the  current  normalized  density  measurement.  To  prepare  for  the  next  iteration  the  co- 
variance  matrix  is  updated  using 
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This  process  of  updating  the  polynomials  of  the  density  estimate  and  then  updating 
the  covariance  matrix  in  preparation  for  the  next  step  is  repeated  at  one  second  intervals 
until  atmospheric  exit.  With  the  density  estimate  now  in  place,  an  estimate  of  the  velocity 
loss  to  occur  in  the  exit  pha.se  due  to  aerodynamic  drag  may  be  computed  by  integrating 
the  drag  equation.  Begin  by  writing  the  drag  equation. 

f  ^  (26) 

d!  2^  '  m  2Mo 
where  Mp  is  the  vehicle  ballistic  coefficient 


u  -  J!L 
CpS‘ 


Replace  dt  with  ^  in  Eq.  (26)  and  substitute  the  expression  given  in  Eq.  (12)  for  p 


and  rearrange  terms  to  obtain 


dV  0.5h{h)  2  6v 

—1  = - ^  (Cj  +  C2Jr  +  C3jr^  + ... +C7J:  )dh. 


Since  x  =  h/h^.  if  h  is  a  function  of  h.  this  expression  can  be  integrated  analytical¬ 
ly  between  any  two  altitudes  to  determine  the  change  in  velocity  due  to  aerodynamic  drag. 
Again,  as  was  done  in  the  ARC  algorithm  and  in  Eq.  (7),  li  is  chosen  to  be  a  constant.  If  h 
is  less  than  250,(XX)  ft,  p  has  a  discontinuity  at  /i,  =  250, 000  ft;  so,  the  integration  must 
be  performed  in  two  steps,  one  from  the  current  altitude  to  250,000  ft  and  a  second  from 
250,0(X)  to  the  exit  altitude.  If  we  change  the  variable  of  integration  from  fi  to  x  we  get 


,  1  0.5p,  Tf*' 

1  = _ ^  {e  LhSlJ(^  + 

V,  M„h  ‘  '  ^ 

X,  D  I  \x, 


0  ft 

c^x  +  ...  +c-jX  )dx 
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(jt-j:,) 


f-l 

LhS2j 


(c ,  +  C2JC  +  c-^x^  + ...  +  c^x^)  dx 


where  x^  =  h^/h^  and  the  current  altitude  is  expressed  x^  =  h/h^.  When  h  is  greater 
than  250,000  ft  this  integration  may  be  carried  out  in  a  single  integration  step. 


1  0.5pj 


1  rM 
r.-<'-'''>Lhs2j 


J 


{c^  +  C2X  +  c-^x^  + ...  +c-jX^)dx 


(30) 


To  integrate  this  expression,  repeatedly  integrate  the  expression  for  density  by  parts 
to  obtain 


-p.e 


-(x-x,) 


hS 


b 


'rhSi  . 

U+  -r-  u'  + 

hS- 

3 

u"  +  ...  + 

hS- 

■'A' 

dx^ . 

(31) 


where 


«  =  c,  +  C2X  +  CjX^  +  c^x^  +  C3/  +  c^x^  +  Cjx^  (32) 

and  the  primes  indicate  a  derivative  with  respect  to  x  so  that 

u'  =  C2  +  2C3X+ 304X^  +  405X^4- +  (33) 

«"  =  203  +  604X  + ISojX^  +  20o^x^  +  30cjX^  (34) 

This  process  is  continued  until  all  six  of  the  required  derivatives  are  formed  using  the  val¬ 
ues  of  0  which  were  most  recently  estimated,  u  and  all  six  of  its  derivatives  are  calculated 
for  the  current  altitude,  and  the  atmospheric  exit  altitude.  Additionally,  u  and  the  six  de¬ 
rivatives  must  be  calculated  at  250,000  ft  altitude  if  the  current  altitude  is  below  250,000 
ft.  These  values  are  inserted  into  Eq.  (31),  which  is  in  turn  inserted  into  Eq.  (29)  or  (30)  as 
required.  The  predicted  velocity  loss  due  to  aerodynamic  drag  is  then  found  by  solving 
Eq.  (29)  or  (30)  for  the  predicted  exit  relative  velocity  and  then  subtracting  the  current  rel¬ 
ative  velocity 
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AV  =  V  -V  . 

rx  r 


If  the  current  altitude  is  below  250,000  ft 


„  _  1  .  .  rhSl-l2 

"  M^h  I  h,  I  h^j 


hSl  ^rhSl-l^  rhSnVul 

-r-u+  -J—  u'  +  ...+  -j—  —r  (36) 

l^e  I  J  dx^ij 


0-5p,  [  -(x-x,)j^rjj52  .rhS2n2 


M^h 


hS2  rhS2-i2  ,  rhSll'^d^u' 

—  u+  —  u  + ...  +  — —  — -T 

K  L  ^  J  L  *  J 


it  "i~* 


or  if  the  current  altitude  is  above  250,000  ft 


1  ,  ®-^Pl  f  "^-"■'''^h^rhS2  ,  rhS2-i2  ,  rhS2nVM 


''rx=  ^  + 

M^h 


hS2  .  rhSZi^  ,  rhS2-i‘ 

-r-M+  -j—  u'  +  ...+  — 

K  IK]  IK] 


Improved  Exit  Phase 

This  improved  exit  phase,  first  published  by  the  Charles  Stark  Draper  Laborato- 

A'X  AA  • 

ry  ,  is  a  simplified  method  of  calculating  the  required  altitude  rate  h  for  the  APC  con¬ 
troller.  It  is  intended  to  replace  the  exit  phase  presented  in  Appendix  B.  To  begin  the 
derivation,  the  velocity  loss  due  to  aerodynamic  drag  is  calculated  starting  with  the  differ¬ 
ential  equation  for  drag; 

T  =  -4-  (38) 

dt  Mjy 

Rearrange  terms  in  Eq.  (38)  to  obtain 

dV  =  -il-d/.  (39) 

Md 

dh 

Replace  dt  with  ~  and  expand  q. 
dh 
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0.5  p 


Moh 


-dh 


(40) 


With  the  assumptions  of  a  constant  altitude  rate  and  an  exponential  atmosphere  of 
known  scale  height  the  above  equation  may  be  integrated  analytically  to  obtain  the  change 
in  velocity  AV  which  will  occur  due  to  aerodynamic  drag.  This  result  is  slightly  different 
from  the  original  APC  exit  phase  derivation  which  uses  this  equation  to  predict  the  veloci¬ 
ty  at  exit  instead  of  computing  the  change  in  velocity  which  will  occur  due  to  drag.  The 
preferred  form  for  the  drag  equation  is 


AV  =  -1/ 


^^des^D  1  'j 

hS^  VgJ- 


(41) 


To  use  the  hybrid  density  estimator  replace  the  expression  for  AV  given  in  Eq.  (41) 
with  the  expression  given  in  Eq.  (10).  Likewise,  to  use  the  polynomial  density  estimator 
replace  the  results  of  Eq.  (41)  with  those  of  Eq.  (35). 

The  desired  velocity  for  a  vehicle  in  a  purely  Keplerian  (no  aerodynamic  forces)  or¬ 
bit  at  the  current  radius  with  the  desired  altitude  rate  h  to  attain  the  targeted  apocenter  ra¬ 
dius  may  be  computed 


des 


^target 

a 

"des 

R{R  + 

R  ^ 

”  target  . 

(42) 


The  first  term  under  the  radical  is  the  velocity  at  pericenter  for  an  elliptical  orbit  with  peri- 
center  radius  R  and  apocenter  at  R^^gc^ 
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target 


A  new  variable  was  introduced 


0.5 


'factor 


f  R  2 

‘  .^target. 

^  R  ^ 


per 


Therefore  may  be  written 


(43) 


(44) 


(45) 


To  avoid  the  square  root  in  Eq.  (45)  a  small  term  is  added  under  the  radical  to  complete  the 
square 


'^des  ^  approximated 


^factor  ^  .2 
^des 


per 


+  0(e^)  = 


J 


•  2 

^pcr  ^factor  ^des 


+  0(e2) 


(47) 


The  corrector  step  to  update  altitude  rate  is  a  single  step  Newton  iteration.  The  dif¬ 
ference  between  the  current  inertial  velocity  minus  the  velocity  loss  expected  from  aero¬ 
dynamic  drag  and  the  desired  velocity  computed  above  is  called  the  velocity  miss  or 

V  ■ 
miss 

=  (Vi-AV)  - '  CS) 
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ilV 

miss  * 

The  negative  of  is  then  divided  by  — ^ —  to  produce  an  update  for  /i(jes  • 

(Jhdes 


^des  update  ~  ^des  ■*" 


( y-^^factor^des)  " 

|j5^  ^factor  ^des 


Equilibrium  Glide  to  Exit  Phase  Transition  Velocity 

To  minimize  total  AV  required  to  transition  from  the  intermediate  orbit  to  the  desired 
orbit  it  is  sufficient  to  minimize  the  exit  flight  path  angle  provided  the  vehicle  exits  in  the 
desired  orbit  plane  and  the  apocenter  of  the  intermediate  orbit  equals  the  desired  apocenter. 
This  approach  will  maximize  the  pericenter  of  the  post-aero  braking  orbit.  If  the  controller 
is  able  to  prof)erly  target  the  apocenter  altitude,  then  minimizing  will  produce  a  maxi¬ 
mum  exit  velocity,  a  maximum  pericenter  for  the  intermediate  orbit  and  a  minimum  AV. 
Fig.  1 1  shows  how  selecting  a  higher  transition  velocity  for  the  APC  controller  to  switch  to 
the  exit  phase  control  algorithm  will  tend  to  minimize  y^  and  the  A  V  required  to  attain  the 
desired  final  orbit  provided  the  vehicle  can  properly  target  the  desired  apocenter.  When  the 
transition  velocity  is  increased  the  predictor/corrector  step  will  calculate  a  lower  li  to  target 
the  desired  apocenter.  The  drawback  to  minimizing  y^  by  increa'^ing  the  transition  velocity 
and  using  a  shallower  flight  path  for  the  exit  phase  is  that  by  doing  so  the  exit  phase  will  be 
flown  using  a  higher  percentage  of  the  available  lift  to  follow  the  desired  trajectory.  In  the 
limit  the  minimum  A  V  path  flies  the  entire  exit  trajectory  with  a  bank  angle  of  1 80° .  When 
the  transition  velocity  becomes  too  great  the  vehicle  can  no  longer  maintain  the  required 
shallow  flight  path,  even  in  a  relatively  smooth  atmosphere,  and  may  overshoot  the  desired 
apocenter  altitude,  as  seen  in  Fig.  12.  Following  this  type  of  shallow  flight  path  angle  tra¬ 
jectory  severely  limits  the  robustness  to  density  dispersions.  If  in  the  initial  phases  of  the 
exit  phase,  the  control  system  calculates  a  shallow  exit  trajectory,  one  whi  'h  requires  al- 
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Fig.  11  Exit  Flight  Path  Angle  and  AV  Required  vs.  Transition  Velocity 


most  full  lift  down  to  maintain,  any  decrease  in  atmospheric  density  from  that  modeled  in 
the  predictor  step  will  result  in  less  control  authority  and  an  inability  to  fly  the  shallow  tra¬ 
jectory,  less  velocity  loss  than  predicted  resulting  in  a  faster  exit  speed  than  desired  and  a 
post-aerobraking  apocenter  higher  than  desired.  An  increase  in  AV  results.  On  the  other 
hand,  transitioning  to  the  exit  phase  at  a  velocity  which  is  too  slow  guarantees  an  increase 
in  A  V  by  requiring  a  steep  h  to  target  apocenter  which  produces  large  exit  flight  path  an¬ 
gles.  The  best  trajectory  is  one  which  strikes  a  desirable  balance  between  minimizing  A  V 
while  retaining  enough  control  to  be  robust  under  the  influence  of  off-nominal  density  vari¬ 
ations.  It  would  seem  to  be  a  simple  matter  to  pick  a  transition  velocity  which  produces  the 
desired  balance,  but  the  “correct”  transition  velocity  varies  with  the  state  of  the  atmosphere, 
the  initial  conditions,  and  the  vehicle  configuration. 
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Fig.  12  Apocenter  and  Pericenter  Altitudes  vs.  Transition  Velocity 


When  selecting  the  appropriate  transition  velocity,  the  important  criteria  are  the  drag 
coefficient  of  the  vehicle  and  the  atmosphere  yet  to  be  traversed.  These  two  parameters 
and  h  define  the  velocity  loss  that  will  occur.  After  considerable  testing  a  desired  altitude 
rate,  h(.,  of  450  ft/sec  was  found  to  yield  a  good  trade  off  between  minimizing  AV  and 
producing  robustness.  The  simulations  depicted  in  Fig.  1 1  and  Fig.  12  require  a  transition 
velocity  of  14,922  ft/sec  to  produce  an  exit  phase  altitude  rate  of  450  ft/sec.  As  seen  in 
Fig.  1 1  this  transition  velocity,  and  hence  this  altitude  rate  are  removed  somewhat  from  the 
region  where  exit  flight  path  angle  and  AV  increase  dramatically.  Yet  this  altitude  rate  was 
still  steep  enough  to  provide  a  measure  of  robustness  against  density  variations.  Armed 
with  this  choice  for  altitude  rate,  a  better  way  to  calculate  transition  velocity  may  be  for¬ 
mulated. 


38 


The  desired  velocity  ''des  for  a  vehicle  in  a  purely  Keplerian  (no  aerodynamic 
forces)  orbit  at  the  current  radius  with  the  desired  altitude  rate  h  to  attain  the  targeted  apo- 
center  radius  was  computed  in  Eq.  (47).  By  adding  the  velocity  loss  expected  due  to  aero¬ 
dynamic  drag  the  current  velocity  required  to  target  the  desired  apocenter  altitude  by 
following  a  450  ft/sec  path  may  be  computed.  The  chosen  method  of  density  estimation 
may  be  used  to  compute  the  velocity  loss  due  to  aerodynamic  drag,  by  inserting  the  de¬ 
sired  450  ft/sec  altitude  rate  into  the  appropriate  derivation.  Equation  (10)  is  used  if  the 
hybrid  density  estimator  is  the  selected  method  of  density  estimation,  whereas,  Eq.  (35)  is 
used  for  the  polynomial  density  estimator  or  Eq.  (41)  for  the  simple  estimate  of  a  constant 
scale  height  exponential  atmosphere.  One  additional  term  is  added  to  allow  for  the  veloc¬ 
ity  loss  between  initiation  of  the  exit  phase  and  achievement  of  the  desired  altitude  rate. 
The  appropriate  velocity  to  transition  from  the  equilibrium  glide  phase  to  the  exit  phase 
may  now  be  expressed 


'',ng;= 


(50) 


5/  in  this  equation  is  the  time  required  from  initiation  of  the  exit  phase  until  the  de¬ 
sired  altitude  rate  is  attained.  The  vehicle  modeled  in  this  study  has  a  limit  of  5°/sec^  on 
roll  acceleration  and  20°/sec  on  roll  rate.  A  value  of  20  seconds  was  selected  for  5r  be¬ 
cause  with  these  current  limits  on  roll  rate  and  roll  acceleration  the  vehicle  requires  thir¬ 
teen  seconds  to  perform  a  180'’  rest  to  rest  maneuver.  After  rolling  to  the  lift  up 
configuration  there  is  still  an  additional  delay  of  five  to  ten  seconds  before  the  vehicle’s  al¬ 
titude  rate  matches  the  desired  value.  With  5/  set  to  20  seconds  the  transition  velocity  cal¬ 
culation  performed  extremely  well.  The  methodology  for  employing  a  variable  transition 
velocity  is  to  compute  using  eq.  (50).  When  the  inertial  velocity  decreases  below 
the  calculated  the  controller  initiates  the  exit  phase. 
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CHAPTER  III 

LYAPUNOV  CONTROLLERS 

Lee  and  Grantham  present  a  Lyapunov  Steepest  Descent  controller  which  they 
claim  is  robust  to  atmospheric  perturbations.  Their  controller  is  for  a  vehicle  which  mod¬ 
ulates  angle  of  attack  while  the  MRSR  vehicle  under  study  flies  at  a  constant  angle  of  at¬ 
tack  and  varies  the  bank  angle  to  control  the  trajectory.  A  similar  controller  is  developed 
to  control  the  MRSR  vehicle.  A  desired  target  state  is  defined  for  the  vehicle  at  atmospher¬ 
ic  exit  which  will  minimize  the  AV  required  to  transition  to  the  desired  final  orbit.  A  pos¬ 
itive  definite  Lyapunov  function  is  defined  such  that  the  vehicle’s  state  is  at  the  target  when 
the  Lyapunov  function  is  zero.  The  control  variable  is  then  selected  so  that  the  Lyapunov 
function  is  driven,  in  a  steepest  descent  fashion,  toward  the  origin.  When  this  method 
failed  to  be  as  robust  as  hoped,  a  new  Lyapunov  Tracking  controller  was  developed. 

The  Lyapunov  Tracking  controller  permits  the  introduction  of  a  preferred  path  lead¬ 
ing  the  vehicle  to  an  exit  state  which  gives  an  acceptable  A  V  to  transition  to  the  desired  fi¬ 
nal  orbit.  In  the  particular  case  studied  here,  the  preferred  path  is  recomputed  for  each 
trajectory  based  on  accelerometer  data  fed  back  to  the  controller  early  in  the  flight  and  a 
“best  guess”  of  the  density  function  for  the  remainder  of  the  trajectory.  Again,  a  positive 
definite  Lyapunov  function  is  defined  such  that,  if  the  vehicle  is  on  the  preferred  path,  the 
Lyapunov  function  is  zero.  The  control  variable  is  again  selected  in  a  “Lyapunov  Opti¬ 
mal”  fashion  to  drive  the  Lyapunov  function  toward  the  origin  as  quickly  as  possible.  A 
gain  scheduling  scheme  defines  an  optimal  descent  function  for  each  phase  of  the  trajecto¬ 
ry,  Finally,  because  of  high  trajectory  loads  generated  by  this  control  scheme  and  difficul¬ 
ty  in  acquiring  the  desired  path,  this  Lyapunov  tracking  controller  was  employed  only 
during  the  exit  phase  following  the  modified  equilibrium  glide^’  phase  of  the  MFC  con- 
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troller  presented  in  Chapter  2.  The  equilibrium  glide  phase  was  developed  to  minimize 
trajectory  loads  and  is  very  good  at  doing  just  that.  With  the  modifications  suggested  in 
Chapter  2  the  equilibrium  glide  phase  control  algorithm  is  very  robust  to  perturbations  in 
atmospheric  density.  The  transition  velocity  from  the  equilibrium  glide  phase  to  this  exit 
phase  is  chosen  using  the  methods  presented  in  Chapter  2  so  the  trajectory  is  at  the  base  of 
the  preferred  path  when  transition  occurs.  The  Lyapunov  Tracking  exit  phase  then  follows 
the  computed  path  to  exit  the  atmosphere  with  exit  state  very  near  the  minimum  AV  exit 
state. 


Lyapunov  Steepest  Descent  Controller 


Equations  of  Motion 

Derivation  of  the  Lyapunov  Steepest  Descent  control  algorithm  begins  with  the 
equations  of  motion  for  planar  flight 


dr  _  dh 
dt  ~ 


Vsiny 


(51) 


dv  ^ 


(52) 


dt 


-Clpsv] 

2mV 


cosO  -  ( 


-  - )  cosy 


(53) 


Eq.  (51)  is  simply  the  radial  velocity  in  terms  of  the  inertial  velocity  and  flight  path 
angle.  Eq.  (52)  gives  the  time  rate  of  change  in  velocity  composed  of  two  parts:  1 )  the  ve¬ 
locity  loss  due  to  aerodynamic  drag  and  2)  the  change  in  velocity  due  to  gravitational  ac¬ 
celeration,  often  referred  to  as  the  inertial  component.  Similarly,  Eq.  (53)  is  the  time  rate 
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of  change  in  the  flight  path  angle,  also  composed  of  two  parts;  1)  the  change  in  flight  path 
angle  due  to  the  component  of  aerodynamic  lift  in  the  vertical  plane  and  2)  the  change  in 
flight  path  angle  due  to  gravitational  acceleration  (the  inertial  component).  The  control 
variable  d>,  the  bank  angle,  determines  the  amount  of  lift  exerted  in  the  vertical  plane  to 
bend  the  trajectory  and  change  the  flight  path  angle. 

Nondimensional  State  Variables 

Dimensionless  state  variables  are  introduced; 


^1 

h/h^ 

X2 

= 

.•"3 

y 

(54) 


along  with  a  dimensionless  time  variable  t 


T  =  (///i  ) 


The  equations  of  motion  may  now  be  written; 


JC|  =  JC2SinjC3 


X2  =  -Bax2  - 


(c-  1  +x,) 


2smx3 


Acsxi 


X-i  = 


-COS<I)  + 


cosjr. 


c  -  l+x. 


ic-l+x^)x2 


(55) 


(56) 

(57) 


(58) 


wherea  =  p/p^  =  exp  [(-(/« -/iq)  )/hS] ,  A  =  (p^S/i^C^)  /  (2ni) , 

D  =  (PqS/i^C^)  /  (2m)  and  c  =  R/h^.  It  has  been  found  that  a  good  approximation  is 
to  assume  that  the  relative  and  inertial  velocity  differ  by  a  constant,  so  that 
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V-&V  (59) 

and  similarly 

X2  =  JC2-5jr2-  (60) 

For  this  controller  it  is  more  convenient  to  replace  cosd*  in  Eq.  (58)  with  the  control 
variable  u  where  u  =  cosO  and  is  thus  bounded  between  ±1.  Eq.  (58)  is  therefore  re¬ 
placed  with 

Aax^u  cosJtj  p  ^ 

^3  ic-l+x^)x2 

Target  State 

The  minimum  A  V  aerobraking  maneuver  is  one  which  exits  the  atmosphere  on  a  tra¬ 
jectory  with  the  correct  apocenter  altitude  and  a  maximum  vacuum  pericenter  altitude. 
This  goal  is  attained  by  exiting  the  atmosphere  with  the  minimum  possible  flight  path  an¬ 
gle  and  the  correct  velocity  to  attain  the  desired  apocenter.  The  goal,  therefore,  is  to  guide 
the  vehicle  along  an  aerobraking  trajectory  which  reaches  the  atmospheric  interface  alti¬ 
tude  with  the  correct  velocity  to  attain  the  desired  apocenter  altitude  while  maintaining  a 
minimum  positive  flight  path  angle  at  exit.  The  flight  path  angle  must  remain  positive  for 
the  vehicle  to  exit  the  atmosphere.  This  design  objective  is  established  by  setting  the  tar¬ 
geted  flight  path  angle  at  atmospheric  exit  to  zero  and  establishing  a  target  exit  velocity. 

^  n 

The  target  state  may  be  pre.sented  in  non-dimensional  form  as 


The  target  exit  velocity,  and  hence  ^2  *«ay  be  derived  assuming  a  Keplerian  orbit 
from  atmospheric  exit  to  apocenter.  This  desired  exit  velocity  is  a  function  of  the  exit 


flight  path  angle,  and  several  constants  for  the  problem  including  the  atmospheric  interface 
radius  (R)  and  target  apocenter  radius  (r^)  ■  The  desired  exit  velocity  is 


Descent  Function 

A  function  is  a  descent  function  if,  and  only  if,  it  is  a  positive  definite  differentiable 
function.  That  is:^^ 


W  (jc)  >0  for  all  jc 

(64) 

0 

II 

(65) 

(66) 

Any  candidate  Lyapunov  function  may  be  chosen  as  the  descent  function  W  [jr  (/)  ] . 
However,  the  most  logical  choice,  and  the  one  recommended  by  Lee  and  Grantham^,  is  a 
weighted  quadratic  measure  of  distance  to  the  target.  This  function  is  expressed 

W{x)  =  -  1  -^2  “^2  ^3]  0  10  -*^2  “-^2 

/?12  Op 22  x-^ 

where  the  constant  weighting  terms  are  chosen  to  define  a  preferred  direction  toward 
the  target  in  the  JCj  -  state  space.  The  preferred  direction  for  the  states  is  presented  in 
Fig.  13.  An  ellipsoid  is  chosen,  oriented  so  that,  while  the  vehicle  is  deep  in  the  atmo¬ 
sphere,  the  preferred  direction  (opposite  the  descent  function  gradient)  in  the  x^  -  x-^  state 
space  gives  positive  lift  to  climb  out  of  the  atmosphere,  but  as  the  vehicle  approaches  at¬ 
mospheric  exit  the  preferred  direction  uses  negative  lift  to  minimize  the  exit  flight  path  an¬ 
gle.  The  weights  must  be  scaled  so  the  velocity  reaches  the  target  velocity  as  the  vehicle 


44 


(Adapted  from  Reference  7) 


reaches  the  atmospheric  interlace  altitude.  For  the  elliptical  descent  function  shown  in 
Fig.  13,  may  be  calculated  as  follows  . 


2  7  2  2 

Pjj  =  a  sin  ^  +  b  cos  (|) 

(68) 

2  2 

Pl2  -  sin<|>cos4> (a  -b  ) 

(69) 

2  2  2  2 

P22  =  b  sin  ({)  +  «  cos  ij) 

(70) 

The  angle  between  the  gradient  of  the  descent  function  and  the  state  space  velocity 
vector  f(x,  u)  is  expressed 


H{x,u) 


OW(x))/(dx)  f{x.u) 

luawcxll/ox)  II  I|/(jr,«)|| 


(71) 


where  P  is  the  angle  between  the  gradient  of  the  descent  function  and  f{x,  u) . 
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Lyapunov  Steepest  Descent  Optimal  Control 

For  the  control  to  be  Lyapunov  steepest  descent  optimal,  u*  (x)  must  make^ 

H[x,u*(x)]  ^H[x,u(x)]  (72) 

for  all  M  6  (/  where  U  is  the  allowable  set  of  controls  bounded  by  ±1.  Furthermore,  for 
the  control  to  be  Lyapunov  steepest  descent  optimal,  f[x,  (j:)  1  *0.  If  it  were  possible 
to  make  H  [x,  u*  (x)  ]  <0  everywhere,  then  global  stability  with  respect  to  the  target 
could  be  guaranteed.  Even  if  H  [x,  u*  {x)  ]  <  0,  asymptotic  stability  with  respect  to  the 
target  could  be  guaranteed.  Unfortunately,  with  u  bounded  between  ±1  neither  of  these  is 
always  possible.  Even  so,  u*  (jc)  tries  to  move  the  system  state  variables  as  nearly  oppo¬ 
site  the  gradient  of  the  descent  function  as  possible,  given  the  dynamics  of  the  system  and 
the  limits  on  the  control. 

To  determine  (x)  set  =  0  and  solve  for  «.  If  this  value  of  u  lies  between  ±1 

then  u*  (of)  is  either  this  value  or  ±1,  whichever  minimizes  H.  If  the  value  of  ti  which 

dH  ,47 

solves  ^  =  0  is  not  between  ±1,  then  u*  (x)  is  selected  from  ±1  to  minimize  . 
au 

Performance  Results 

The  Lyapunov  Steepest  Descent  feedback  control  algorithm  will  guide  the  vehicle 
to  very  near  the  minimum  AV  exit  state  provided  the  weights,  and  hence,  W(x)  is 
properly  selected.  Unfortunately,  those  weights  must  be  readjusted  to  attain  accept¬ 
able  performance  for  each  perturbed  entry  condition,  vehicle  lift  and  drag  perturbation,  or 
atmospheric  density  perturbation.  No  acceptable  method,  other  than  a  manual  search,  was 
found  to  determine  the  appropriate  weighting  for  each  perturbed  run.  Clearly,  this  lack  of 
asymptotic  stability  is  not  compatible  with  the  objectives  of  this  research. 
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An  appropriate  descent  function  was  found  for  this  controller  to  perform  acceptably 
for  the  nominal  case  of  a  vehicle  targeting  a  270  nm  circular  orbit  after  entering  the  Mar¬ 
tian  standard  atmosphere  with  6.0  km/sec  velocity  relative  to  the  planet,  -12°  entry  flight 
path  angle,  and  a  lift  to  drag  ratio  of  I  and  for  the  same  vehicle  and  entry  conditions  en¬ 
countering  a  high  or  low  density  Martian  atmosphere.  The  ellipse  which  determines  the 
descent  function  was  chosen  to  have  a  semimajor  axis  of  1.65,  a  semiminor  axis  of  0.41, 
and  a  rotation  angle  <])  of  4.2°,  measured  as  shown  in  Fig.  13.  A  few  perturbations  from 
the  nominal  case  are  then  simulated  and  the  somewhat  disastrous  results  are  presented  in 
Table  1  along  with  the  optimal  results  for  the  same  perturbations  generated  using  the 
method  of  Appendix  A. 

This  controller  was  not  as  robust  as  hoped,  given  the  density,  navigation  or  vehicle 
perturbations  expected  for  the  Martian  aerobraking  problem.  The  controller  may  be  fine 
tuned  for  one  rate  of  energy  depletion,  but  if  anything  alters  the  rate  of  energy  loss  the  con¬ 
troller  must  be  readjusted,  by  altering  the  relative  weights  between  the  states,  to  bring  the 
velocity  to  the  targeted  velocity  just  as  the  vehicle  passes  through  the  atmospheric  inter¬ 
face  altitude.  A  steeper  entry  flight  path  angle  will  thrust  the  vehicle  deeper  into  the  atmo¬ 
sphere,  thereby  increasing  the  rate  of  energy  loss.  Likewise  an  atmosphere  which  is  more 
dense  than  expected,  or  a  drag  coefficient  higher  than  expected,  will  cause  the  vehicle  to 
lose  energy  at  a  higher  rate  than  planned,  resulting  in  exit  conditions  which  are  too  slow 
and  an  apocenter  altitude  which  is  too  low.  In  the  worst  instances  the  vehicle  fails  to  exit 
the  atmosphere  at  all.  Similarly,  a  shallower  entry  flight  path  angle,  less  dense  atmosphere 
or  lower  drag  coefficient  will  result  in  less  velocity  loss  than  needed  and  apocenter  alti¬ 
tudes  higher  than  desired.  To  reduce  the  sensitivities  to  perturbations  which  change  the 
rate  of  energy  loss  the  Lyapunov  controller  was  reformulated  as  a  tracking  controller  de¬ 
signed  to  follow  a  chosen  path  to  atmospheric  exit. 


Table  1  Lyapunov  Steepest  Descent  Controller  Performance  Results 
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Lyapunov  Tracking  Controller 

In  an  effort  to  gain  acceptable  robustness  for  this  controller  the  methodology  was 
changed  from  a  steepest  descent  controller  which  targets  the  optimal  terminal  state  to  a 
steepest  descent  controller  which  targets  a  preferred  path.  That  path  then  is  selected  to 
lead  the  vehicle  to  a  desirable  exit  state  with  enough  robustness  to  prevent  minor  density 
upsets  from  being  catastrophic. 

The  Preferred  Path 

Derivation  of  this  controller  begins  with  definition  of  the  preferred  path.  As  with  the 
predictor  corrector  algorithms  a  constant  altitude  rate  path  leading  to  the  desired  atmo¬ 
spheric  exit  state  was  selected.  The  difference  between  this  Lyapunov  Tracking  Controller 
(LTC)  and  the  predictor  correctors  is  in  how  the  controller  computes  the  constant  altitude 
rate  path.  The  predictor  corrector  algorithms  use  various  methods  to  select  a  constant  alti¬ 
tude  rate  which  will  give  the  desired  apocenter  altitude  and  then  use  altitude  rate  error  to 
select  the  appropriate  bank  angle.  The  LTC,  on  the  other  hand,  assumes  that  it  is  desirable 
to  always  fly  the  same  altitude  rate  to  atmospheric  exit  and  arrive  there  with  the  appropri¬ 
ate  velocity  to  achieve  the  proper  apocenter  altitude.  The  LTC  then  selects  the  in  plane 
portion  of  lift  to  approach  the  path  in  a  steepest  descent  fashion.  The  chosen  altitude  rate 
is  selected  to  produce  the  desired  trade-off  between  robustness  to  density  perturbations 
and  minimum  AV. 

A  constant  altitude  rate  of  450  ft/sec  was  again  selected  (as  on  page  37)  to  define  the 
desired  path  leading  to  atmospheric  exit  with  the  appropriate  velocity  to  target  the  desired 
apocenter  altitude.  This  altitude  rate  produces  trajectories  which  require  within  20  to  30 
ft/sec  of  the  minimum  A  V  values  for  the  various  expected  perturbations  without  short  pe¬ 
riod  density  upsets,  yet  is  still  robust  to  density  variations  of  ±50%  over  small  altitude  in- 
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tervals.  The  equations  used  to  derive  the  i  .iproved  exit  phase  for  the  Mars  Predictor 
Corrector  are  employed  again  here  to  define  this  path.  The  velocity  required  at  a  given  al¬ 
titude  flying  a  specified  altitude  rate  assuming  a  Keplerian  orbit  (no  aerodynamic  drag  ef¬ 
fects)  was  given  in  Eq.  (47)  but  is  repeated  here  for  completeness. 


V.  =  V 
des  per 


1  + 


^factor  ^^2 

des 


per 


P 


=  V. 


per  ^factor  ^des 


(73) 


The  velocity  loss  expected  due  to  aerodynamic  drag  must  be  added  to  this  velocity  to 
determine  the  current  •  clocity  for  the  desired  path.  Note  that  this  desired  velocity  is  a 
function  of  the  dynamics  of  the  Martian  orbit,  the  current  altitude,  the  selected  altitude  rate 
(450  ft/sec)  and  the  expected  velocity  loss  (which  is  a  function  of  the  expected  atmospher¬ 
ic  density  function  and  the  vehicle  coefficient  of  drag).  The  velocity  loss  expected  due  to 
aerodynamic  drag  is  calculated  assuming  the  450  ft/sec  altitude  rate  path  will  be  flown  us¬ 
ing  Eq.  (10)  of  the  hybrid  density  estimator,  or  with  Eq.  (35)  of  the  polynomial  density  es¬ 
timator,  or  with  Eq.  (41)  using  the  simplification  of  a  constant  scale  height  exponential 
atmosphere.  The  desired  current  velocity  defining  the  preferred  path  is 

(74) 

This  velocity  may  be  converted  to  non-dimensional  form 

=  C'/(7|i7r)  (75) 


0=  + 


The  desired  flight  path  angle  JC3  is  computed 

X3  =  asin  (/jjes/VO  .  (76) 

Together,  ^2  ^nd  ^2  define  the  preferred  path  which  will  lead  the  vehicle  along  a  ro¬ 
bust  corridor  to  a  desirable  exit  state.  Now,  a  Lyapunov  function  must  be  formulated  and  a 
control  found  which  will  drive  the  vehicle  onto  and  then  down  the  chosen  path. 
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The  Lyapunov  Descent  Function 


The  selected  positive  definite  Lyapunov  function  is 


Pll  Pi2 

H-h 

P\2  Pll 

(77) 


This  function  is  analogous  to  distance  from  the  target  path  and  is  zero  whenever  the  vehi¬ 
cle  is  on  the  target  path  and  positive  otherwise.  Again,  the  values  are  chosen  to  form 
an  ellipsoid,  the  negative  gradient  of  which  defines  the  preferred  approach  to  the  target 
path.  This  ellipsoid  is  shown  in  Fig.  14.  The  values  are  computed  from  the  semima- 


Fig.  14  Lyapunov  Tracking  Controller  X2  -  ^^3  Descent  Function 
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jor  axis  a,  the  semitninor  axis  b,  and  the  angle  of  rotation  (j)  of  the  ellipse  defining  the  pre¬ 
ferred  gradient^  onto  the  chosen  path  as  in  Eqs.  (68)  through  (70).  They  are  repeated  here 
for  completeness 


Pjj  =  a^sin^4i +  6^cos^(|)  (78) 

Pj2  =  sin<l)cos(l)(n^-^^)  (79) 

P22  ~  ^^sin^ij) +  a^cos^(j)  (80) 


Again,  the  angle  between  the  gradient  of  the  descent  function  and  f(x,  u)  is  ex¬ 
pressed 


{dW{x))/idx)  f{x,u) 
||OlV(jt))/Ox)||||/(x.M)|| 


(81) 


where  P  is  the  angle  between  the  gradient  of  the  descent  function  and  the  state  space  ve¬ 
locity  vector  f(x,u) . 


Selection  of  the  Control 

As  was  done  for  the  Lyapunov  Steepest  Descent  Controller,  a  control  u*  (x)  is 
sought  which  will  move  the  system  state  variables  as  nearly  opposite  the  gradient  of  the 
descent  function  as  possible,  given  the  dynamics  of  the  system  and  the  limits  on  the  con¬ 
trol. 

3// 

As  before,  to  determine  u*  (x)  set  ^  =0  and  solve  for  u.  Note  however,  that  H  in 

this  discussion  is  not  the  same  function  as  H  in  the  LSDC  discussion.  If  the  value  of  u  lies 

between  ±1  then  u*  (x)  is  either  this  value  or  ±1,  whichever  minimizes  H.  If  the  value  of 
dH 

u  which  solves  3-  =  0  is  not  between  ±1  then  u  (x)  is  selected  from  ±1  to  minimize  H. 
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Performance  Results 

An  acceptable  descent  function  for  this  algorithm  was  not  found.  Fig.  15  shows  the 
X2  and  X2  time  histories  while  Fig.  16  shows  the  x^  and  ^3  time  histories  for  a  trajectory 
guided  by  this  LTC.  The  elliptical  descent  function  was  chosen  to  have  a  semimajor  axis 
of  40,  a  semiminor  axis  of  1,  and  a  rotation  angle  (j>  of  19°,  measured  as  shown  in  Fig.  14. 
This  descent  function  produced  an  apocenter  altitude  following  the  aerobraking  maneuver 
of  356  nautical  miles  which  was  as  close  to  the  270  nm  target  as  possible  while  exiting 
with  a  flight  path  angle  near  the  optimal.  But  as  Fig.  15  shows,  once  the  velocity  neared 
the  target  path  it  failed  to  close  in  and  make  the  final  correction  necessary.  A  second  de¬ 
scent  function  was  found  which  would  guide  the  vehicle  closer  to  the  target  apocenter  alti¬ 
tude.  This  descent  function  used  an  elliptical  function  with  the  same  semi-major  axis  of 
40,  semi-minor  axis  of  1,  but  the  rotation  angle  was  changed  to  55°.  As  Fig.  17  and 
Fig.  18  plainly  show,  the  desired  apocenter  altitude  was  not  achieved  by  following  the  de¬ 
sired  path  to  exit  but  rather  by  reducing  the  velocity  (^2)  more  than  desired  and  then 
climbing  with  a  steeper  flight  path  angle  (jr3)  than  preferred.  Almost  by  accident,  the  de¬ 
sired  apocenter  altitude  was  attained.  No  fixed  configuration  for  the  elliptical  descent 
function  was  found  which  would  allow  the  algorithm  to  acquire  the  target  path  and  follow 
it  to  an  acceptable  exit  state. 

Intuitively,  it  is  easy  to  see  that,  if  the  velocity  is  slower  than  V  and  the  flight  path 
angle  is  less  than  ^3,  the  logical  choice  is  to  use  positive  lift  to  get  closer  to  the  path.  Like¬ 
wise,  if  V  is  greater  than  ^  and  the  flight  path  angle  is  greater  than  Jf3  a  lift  down 
orientation  is  required  to  approach  the  path.  The  ambiguous  areas  are  in  the  other  two 
quadrants  where  either  the  velocity  is  too  fast  yet  the  flight  path  angle  is  too  shallow,  or 
where  the  velocity  is  too  slow  but  the  flight  path  angle  is  greater  than  desired.  It  is  desir¬ 
able  to  define  a  line  passing  through  the  target  state  at  each  instant  in  tinie  separating  the 
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selected  was  the  line  drawn  in  the  X2  -  plane  through  (Jt2,  JP3)  by  solving  for  the  de¬ 
sired  X2  and  x^  states  using  a  slightly  different  altitude  rate.  Though  the  chosen  450  ft/sec 
altitude  rate  works  very  well,  it  should  still  be  acceptable  to  fly  a  425  or  a  475  ft/sec  alti¬ 
tude  rate  path  which  would  lead  the  vehicle  to  the  desired  apocenter  orbit. 


To  define  this  switching  line  the  angle  of  rotation,  ({>,  of  the  descent  function  shown 
in  Fig.  14  was  varied  during  the  trajectory  such  that 


<t>  =  atan 


-  <£^3- 

dt2 
.  . 


(82) 


In  effect  this  defines  a  switching  line  formed  by  linearizing  about  the  current  target  state 
and  varying  altitude  rate.  Though 


dx 


(83) 


because  of  the  missing  ^  ^  component.  But,  since  ^  is  small  compared  to  the  ele- 

ments  of  ^  and  (|>  varies  slowly,  increasing  monotomically  from  about  15°  to  about  75° 
ax 

during  the  exit  phase,  this  component  was  assumed  to  be  insignificant.  The  commanded 
bank  angle  was  still  determined  as  before  by  selecting  the  value  of  u  which  minimizes 
H(x,u)  with  H(x,u)  defined  as  before  in  Eq.(81). 


Though  this  method  of  varying  the  weighting  matrix  showed  improvement,  the  algo¬ 
rithm  still  had  problems  acquiring  the  target  path.  The  vehicle  still  used  lift  down  too  ear¬ 
ly  and  plunged  deeply  into  the  atmosphere,  creating  extremely  high  vehicle  accelerations 
and  heat  rates  in  the  process.  To  cure  this  problem  the  Lyapunov  Tracking  Algorithm 
(LTA)  developed  here  was  incorporated  as  an  exit  phase  following  the  equilibrium  glide 
phase  of  the  MFC  algorithm  presented  in  Chapter  2. 
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Lyapunov  IVacking  Controller  Exit  Phase 

The  equilibrium  glide  phase  was  developed  to  guide  the  vehicle  into  the  atmosphere 
and  hold  the  vehicle  in  equilibrium  until  the  velocity  has  been  appropriately  reduced.  It 
was  designed  to  perform  this  task  while  keeping  the  maximum  trajectory  loads  and  peak 
heat  rates  to  a  minimum.  It  performs  this  task  very  well.  On  the  other  hand  the  LTC  just 
described  performs  well  in  holding  the  desired  path  to  exit  if  somehow  it  could  be  started 
near  that  path.  The  marriage  of  the  LTA  as  an  exit  phase  with  the  equilibrium  glide  phase 
was  implemented  next.  The  method  of  computing  transition  velocity  from  the  equilibrium 
glide  phase  to  the  exit  phase  presented  in  Chapter  2  placed  the  vehicle  very  near  the  de¬ 
sired  path.  This  combination  of  equilibrium  glide  phase  and  LTA  proved  to  be  the  best 
controller  examined. 

The  two  density  estimation  techniques  presented  in  Chapter  2  were  also  tested.  The 
complete  control  algorithm  with  the  hybrid  density  estimator  included  will  be  referred  to 
as  the  Hybrid  Lyapunov  Tracking  Controller  (HLTC)  while  this  controller  with  the  poly¬ 
nomial  density  estimator  will  be  called  the  Lyapunov  Tracking  Controller  (LTC). 
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CHAPTER  IV 

CONTROLLER  SENSITIVITY  ANALYSIS 

The  performance  of  the  MPC  and  MHPC  control  algorithms  developed  in  Chapter  11 
and  of  the  LTC  and  LHTC  control  algorithms  developed  Chapter  111  was  determined  along 
with  the  performance  of  the  APC  and  Energy  Controllers  presented  in  Appendix  B  and 
Appendix  C  respectively.  The  algorithms  were  tested  using  the  six  degree  of  freedom 
computer  simulation  based  on  the  Program  to  Optimize  Simulated  Trajectories^^  (POST), 
which  uses  a  fourth  order  Runge-Kutta  numerical  integration  scheme  to  continuously  inte¬ 
grate  both  the  force  and  moment  equations  of  the  vehicle.  The  control  algorithms  were 
tested  to  determine  the  effect  of  large  scale  density  variations  such  as  those  caused  by  the 
seasonal  sublimation  and  condensation  of  the  Martian  atmosphere  or  by  a  global  dust 
storm.  They  were  also  tested  to  determine  the  effect  of  short  period  atmospheric  varia¬ 
tions  by  injecting  square  wave  density  pulses,  similar  to  those  used  by  Fitzgerald*®’  *  *,  of 
various  magnitudes  and  durations  into  the  density  function  at  various  altitudes.  Entry 
flight  path  angles  were  varied  within  the  current  predicted  error  band"*^.  Perturbations  in 
the  vehicle  lift  and  drag  characteristics  were  also  simulated.  Finally,  combinations  of 
these  perturbations  in  the  atmospheric  density  function,  entry  flight  path  angle  and  vehicle 
lift  and  drag  characteristics  were  simulated  and  the  performance  of  each  controller  was  de¬ 
termined. 

Following  a  brief  description  of  the  vehicle  and  trajectory  simulation  program  used, 
the  data  from  this  test  program  are  presented  graphically  utilizing  three  dimensional  mesh 
plots.  The  primary  thrust  of  this  test  program  was  to  select  the  best  controllers  from  those 
studied.  A  full  performance  evaluation  of  the  selected  controllers,  aimed  at  determining 
the  robustness  limits  of  the  selected  controllers,  is  presented  in  Chapter  V. 
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Vehicle  and  IVajectory  Simulation  Inputs 

The  vehicle  used  in  the  study  is  a  biconic  aeroshell  design  with  a  fifteen  foot  base  di¬ 
ameter  and  a  weight  of  1 1,023  lbs.  The  base  surface  area  of  176. 1  ft^  is  used  as  the  refer¬ 
ence  surface  area.  The  vehicle  is  designed  with  a  five  foot  eg  offset  resulting  in  a  trim 
angle-of-attack  of  27°  which  is  maintained  throughout  the  maneuver  via  a  simple  propor¬ 
tional  feedback  controller.  Control  is  via  bank  maneuvers  which  reorient  the  direction  of 
the  lift  vector.  These  bank  maneuvers  are  commanded  as  body  axis  rolls  with  coordinat¬ 
ing  body  axis  yaw  maneuvers.  The  nominal  lift  coefficient  is  0.68892  while  the  drag  coef¬ 
ficient  is  0.69819,  producing  a  nominal  L/D  of  0.99. 

The  Mars  Global  Reference  Atmosphere  ModeP^  (MARS-GRAM)  was  used  to  pro¬ 
duce  realistic  atmospheres  for  the  study.  Three  different  atmospheres  representing  a  nom¬ 
inal,  a  low  density  and  a  high  density  Martian  atmosphere  were  used  (Fig.  19).  The 
nominal  atmosphere  is  the  COSPARV  Model  Atmosphere  For  Mars'*,  while  the  high  den¬ 
sity  and  low  density  atmospheres  were  derived  using  MARS-GRAM.  The  low  density  at¬ 
mosphere  is  a  MARS-GRAM  simulation  of  the  lowest  density  Martian  atmosphere 
predicted  for  April  10,  1999  assuming  no  dust  storms  and  a  10.7  cm  solar  flux  of  50  (nom¬ 
inal  value  =  150).  The  high  density  atmosphere  represents  the  highest  density  atmosphere 
predicted  on  December  27,  1997,  again  with  no  dust  storms  but  this  time  with  a  10.7  cm 
solar  flux  of  300.  Although  MARS-GRAM  was  incorporated  as  a  subroutine  to  POST 
which  can  be  called  to  generate  atmospheric  data  on  line,  MARS-GRAM  was  not  utilized 
in  this  manner  because  of  the  added  computational  time.  MARS-GRAM  was  used  to  gen¬ 
erate  atmospheric  data  which  were  stored  in  tabular  form.  These  tables  of  atmospheric 
data  were  then  included  in  the  POST  input  namelist. 

In  addition  to  the  large  scale  density  variations  introduced  by  using  the  low,  nominal 
or  high  density  atmosphere  models  described  above,  short  period  variations  in  the  atmo- 
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Fig.  19  Mars  Nominal,  Low  and  High  Density  Atmospheres 


spheric  density  function  were  investigated  by  introducing  square  wave  density  pulses,  sim¬ 
ilar  to  those  used  by  Fitzgerald*®’  *  * .  These  pulses  perturb  the  local  atmosphere  within  a 
10,000  or  20,000  ft  altitude  band  by  multiplying  the  expected  density  by  a  constant  mag¬ 
nitude  density  multiplier.  The  magnitudes  of  the  density  multipliers  used  include  0.5, 
0.75, 1 .25, 1 .5,  1 .75  and  2.0.  The  lower  edge  of  the  density  pulses  were  varied  in  10,000  ft 
steps  from  100,000  to  290,000  ft 

The  atmospheric  interface  altitude  was  selected  to  be  125  km  (410,105  ft)  and  the 
initial  conditions  are  defined  at  this  altitude.  The  entry  velocity  is  6  km/sec  (19,685  ft/sec) 
and  the  nominal  entry  flight  path  angle  is  -12°.  The  targeted  orbit  is  a  270  nm  circular  or¬ 
bit.  In  addition  to  the  atmospheric  perturbations  mentioned  above,  perturbations  were  in¬ 
troduced  in  the  vehicle  lift  and  drag  coefficients  representing  variations  of  ±33%  from  the 
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nominal  LTD  ratio  of  0.99.  The  +33%  L/D  perturbation  was  introduced  by  multiplying  the 
nominal  lift  coefficient  by  1.14,  while  the  drag  coefficient  was  multiplied  by  0.86.  The 
-33%  L/D  perturbation  was  introduced  by  multiplying  the  nominal  lift  coefficient  by  0.8, 
while  the  drag  coefficient  was  multiplied  by  1 .2.  This  method  of  varying  L/D  also  per¬ 
turbed  the  ballistic  coefficient  of  the  vehicle.  Navigation  errors  in  the  form  of  variations  in 
the  entry  flight-path-angle  of  ±0.25°  and  ±0.5°  from  the  nominal  -12°  were  considered. 
The  performance  for  each  perturbed  run  is  presented  as  total  AV  required  to  achieve  the 
desired  final  orbit.  AV  is  a  measure  of  the  controllers  overall  success  in  meeting  the  de¬ 
sired  exit  conditions.  The  AV  was  calculated  assuming  one  burn  at  atmospheric  exit  ori¬ 
ented  along  the  velocity  vector  to  correct  any  apocenter  error,  a  second  at  apocenter  to 
raise  pericenter,  and  a  final  bum  to  correct  any  orbit  plane  error. 

Analytic  Predictor  Corrector  Performance  Pesults 

The  original  APC  controller  presented  in  Appendix  B  did  not  fare  very  well  when 
challenged  with  the  possible  perturbations  used  in  this  study.  Fig.  20,  Fig.  21,  Fig.  22  and 
Fig.  23  present  the  performance  of  this  controller  when  faced  with  these  perturbations. 
Fig.  20,  Fig.  21  and  Fig.  22  show  the  performance  of  a  nominal  vehicle  which  enters  the 
atmosphere  with  an  =  -12°  and  then  encounters  a  square  wave  density  pulse.  AV  re¬ 
quired  to  circularize  is  plotted  on  the  vertical  axis.  Fig.  20  presents  the  results  of  encoun¬ 
tering  square  wave  density  pulses  in  a  nominal  Martian  atmosphere,  while  Fig.  21  presents 
the  results  for  a  low  density  atmosphere  and  Fig.  22  presents  the  results  for  a  high  density 
atmosphere.  In  the  first  diagram  of  each  figure  the  density  pulse  perturbs  a  10,000  ft  alti¬ 
tude  band  while  in  the  second  diagram  the  pulse  affects  a  20,000  ft  band.  Magnitudes  for 
these  pulses  range  from  -50%  to  +100%  in  25%  increments.  The  location  of  the  lower 
edge  of  the  pulse  was  moved  from  100,000  ft  to  290,000  ft  in  10,000  ft  increments.  The 


Fig.  23  Analytic  Predictor  Corrector  Sensitivity  to  Variations  in  Lift  to  Drag 
Ratio  and  Entry  Flight  Path  Angle,  a)  Nominal  Atmosphere;  b)  Low  Density 
Atmosphere;  c)  High  Density  Atmosphere 
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magnitude  of  the  pulse  and  the  altitude  of  the  lower  edge  are  shown  on  the  diagrams.  Fig. 
23  presents  AV  required  to  circularize  as  L/D  and  are  varied.  Fig.  23a  shows  the  results 
of  these  variations  in  a  nominal  atmosphere  while  Fig.  23b  presents  the  same  perturbations 
in  a  low  density  atmosphere  and  Fig.  23c  is  in  a  high  density  atmosphere. 

Fig.  23  shows  that  the  APC  controller  exhibits  a  considerable  sensitivity  to  off  nom¬ 
inal  vehicle  design  and  to  navigation  errors.  The  controller  also  shows  a  marked  decrease 
in  performance  in  the  high  density  atmosphere  with  no  density  steps  when  compared  to  its 
performance  in  the  low  density  and  nominal  atmospheres.  The  AV  required  to  circularize 
following  the  aerobraking  maneuver  in  a  high  density  atmosphere  even  with  no  density 
step  (0%  magnitude  density  step)  is  over  400  ft/sec  while  the  optimal  results  presented  in 
Table  1  show  that  it  should  require  less  AV  to  circularize  after  aerobraking  in  a  high  densi¬ 
ty  atmosphere  (optimally  about  316  ft/sec)  than  in  a  nominal  or  low  density  atmosphere. 
Additionally,  the  variations  in  AV  shown  in  Fig.  23c  are  considerably  worse  than  those  in 
Fig.  23a  or  b.  Part  of  this  sensitivity  comes  from  using  a  specified  transition  velocity  to 
switch  to  the  exit  phase,  ignoring  the  actual  energy  loss  to  occur  during  the  exit  phase. 
The  other  reason  for  this  sensitivity  is  the  rather  simplistic  density  model.  However,  the 
controller  is  less  sensitive  to  density  steps  in  the  high  density  atmosphere  than  in  the  low 
or  nominal  atmosphere.  This  sensitivity  can  again  be  explained  by  the  choice  of  a  speci¬ 
fied  transition  velocity  for  the  switch  from  entry  to  exit  phase. 

The  transition  velocity  selected  for  this  controller  was  14,922  fVsec.  This  transition 
velocity  is  appropriate  for  a  nominal  vehicle  which  enters  the  nominal  atmosphere  with  a 
flight  path  angle  of  -12°.  However,  if  the  initial  flight  path  angle  is  steeper  than  -12°  or 
the  atmosphere  is  more  dense  than  expected,  the  vehicle  will  plunge  into  the  atmosphere 
deeper  than  expected,  and  consequently,  will  have  more  atmosphere  to  traverse  during  the 
exit  phase  and  will  loose  more  energy  to  aerodynamic  drag.  Similarly,  if  the  vehicle’s 
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drag  coefficient  is  higher  than  planned  for,  the  vehicle  will  loose  more  energy  during  exit 
than  planned  for.  Forcing  the  vehicle  to  decelerate  to  a  predetermined  velocity  before  ini¬ 
tiating  the  exit  phase  requires  the  exit  phase  to  be  flown  at  a  steeper  altitude  rate  than  de¬ 
sired  to  target  the  desired  apocenter  altitude  resulting  in  higher  AV  values  and  tends  to 
amplify  the  effect  of  off  nominal  entry  condition  or  drag  coefficient.  The  steeper  exit  path 
flown  by  this  controller  in  the  high  density  atmospheres  however,  is  more  robust  to  density 
variations  as  may  be  seen  in  Fig.  22.  This  result  is  due  to  the  fact  that  a  trajectory  which 
flies  a  steeper  exit  phase  has  reduced  more  of  the  velocity  deep  in  the  atmosphere  and  is 
not  requiring  as  much  velocity  loss  during  the  exit  phase.  Density  variations  which  per¬ 
turb  the  amount  of  velocity  loss  which  actually  occur  during  the  exit  phase  have  less  effect 
when  more  of  the  velocity  is  reduced  deep  in  the  atmosphere.  Later,  the  effect  of  making 
this  transition  velocity  an  adaptive  parameter  will  be  shown. 

The  second  area  of  concern  with  this  controller  is  the  density  estimation  technique. 
The  density  estimator  built  into  this  algorithm  assumes  the  density  function  is  a  fixed  scale 
height  exponential  function.  The  density  derived  onboard  from  accelerometer  measure¬ 
ments  is  filtered  using  a  low  pass  filter  to  remove  high  frequency  noise.  The  result  is  then 
used  to  bias  the  exponential  function  used  to  estimate  density.  This  technique  works  well 
as  long  as  the  density  function  does  not  vary  much  from  a  smooth  exponential  function 
and,  more  critically,  the  scale  height  of  the  atmosphere  is  fairly  constant  and  doesn’t  vary 
much  from  the  assumed  scale  height.  Unfortunately,  the  scale  height  of  the  Martian  atmo¬ 
sphere  does  vary  considerably.  Fig.  19  shows  the  range  in  the  density  function  predicted 
by  MARS-GRAM.  This  figure  presents  altitude  versus  log  density.  The  scale  height  may 
be  determined  by  taking  the  negative  of  the  slope  of  the  density  function  from  this  graph. 
The  scale  height  does  not  vary  considerably  below  250,000  ft,  but  above  250,000  ft  there 
is  considerable  variation.  This  variation  does  not  affect  the  trajectories  flown  in  the  low 
density  atmosphere  very  much  because,  in  a  low  density  atmosphere  above  250,000  ft. 
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aerodynamic  forces  have  negligible  effect  on  the  vehicle.  But  for  those  trajectories  flown 
in  the  high  density  atmosphere  failure  to  properly  model  the  atmosphere  has  considerable 
effect.  The  estimated  density  falls  short  of  the  actual  values  above  250,000  ft.  Conse¬ 
quently,  there  is  more  aerodynamic  drag  than  predicted.  This  typically  manifest  itself  as  a 
final  apocenter  altitude  ten  to  twenty  nautical  miles  lower  than  targeted.  This  problem 
compounded  by  the  steep  altitude  rate  in  the  exit  phase,  brought  on  by  the  constant  transi¬ 
tion  velocity,  lowered  the  ability  of  the  control  system  to  correct  for  density  upsets  which 
occur  after  the  exit  phase  is  initiated  causing  the  sensitivities  shown  in  the  square  wave 
density  pulse  data  of  Fig.  20,  Fig.  21  and  Fig.  22. 

Overall,  the  APC  is  still  a  good  controller.  It  guides  the  vehicle  through  the  aero- 
braking  maneuver  with  minimal  heat,  acceleration,  and  dynamic  pressure  loads,  exiting 
with  an  exit  state  near  optimal  when  the  density  function  encountered  is  near  the  nominal 
value,  when  navigation  is  good  enough  to  allow  precise  control  over  the  entry  state,  and 
when  the  hypersonic  lift  and  drag  characteristics  of  the  vehicle  are  close  to  the  design  val¬ 
ues.  As  part  of  the  modification  of  this  controller  to  meet  the  Martian  requirements  the 
value  of  K-  was  changed  to  4.5  as  recommended  in  Chapter  2.  This  kept  the  vehicle  from 
exiting  the  atmosphere  before  slowing  enough  to  transition  to  the  exit  phase  (skipping  out) 
for  all  of  the  test  cases  examined.  However,  the  APC  just  is  not  quite  robust  enough  to  ad¬ 
equately  handle  the  expected  perturbations  in  the  Martian  aMnosphere,  vehicle  entry  con¬ 
ditions,  or  vehicle  lift  and  drag  variations. 

Energy  Controller  Performance  Results 

Fig.  24,  Fig.  25,  Fig.  26  and  Fig.  27  show  the  performance  of  the  Energy  Controller. 
Again,  the  results  are  pre.sented  in  the  same  format  as  before  with  Fig.  24  showing  results 
for  density  steps  in  a  nominal  atmosphere  while  Fig.  25  shows  the  results  for  a  low  density 
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atmosphere  and  Fig.  26  shows  the  performance  results  in  a  high  density  Martian  atmo¬ 
sphere.  Fig.  27  shows  the  effect  of  varying  Lift  and  Drag  and  entry  flight  path  angle  with 
Fig.  27a  being  in  a  nominal  atmosphere  while  Fig.  27b  show  these  results  in  a  low  density 
atmosphere  and  Fig.  27c  shows  the  results  if  a  high  density  atmosphere  is  encountered. 

The  Energy  Controller  is  substantially  more  robust  than  the  APC  controller  with  re¬ 
spect  to  vehicle  lift  and  drag  perturbations,  and  to  navigation  errors.  Fig.  27a,  b  and  c  all 
show  practically  no  variation  in  AV  required  to  circularize.  Furthermore,  these  results  all 
fall  below  400  ft/sec  to  circularize.  This  insensitivity  may  be  attributed  to  the  fact  that  the 
Energy  Controller  does  not  assume  a  density  function,  though  an  exponential  function  is 
expected;  instead,  it  reties  on  the  current  energy  rate  and  energy  error  to  determine  which 
path  should  be  pursued.  Variations  in  the  vehicle’s  drag  coefficient  simply  changes  the  en¬ 
ergy  rate  and  the  controller  compensates  for  this.  Likewise,  variations  in  the  overall  state 
of  the  atmosphere  (low,  nominal,  or  high  density  atmosphere),  or  variations  in  the  entry 
flight  path  anple  which  force  the  vehicle  deeper  or  shallower  into  the  atmosphere  are  seen 
by  the  controller  as  changes  in  the  energy  rate.  Since  the  controller  seeks  to  make  energy 
rate  approach  zero  as  energy  error  approaches  zero,  variations  of  this  type  are  handled 
well. 


This  method  works  well  as  long  as  the  density  function  is  a  smooth  exponential  but, 
as  the  10,000  and  20,000  ft  density  pulse  diagrams  illustrate,  the  Energy  Controller  shows 
definite  sensitivity  to  density  functions  which  are  not  smooth.  The  large  magnitude 
20,000  ft  duration  density  steps  prove  to  be  more  than  this  controller  can  tolerate.  Fig.  25b 
shows  that  in  a  low  density  atmosphere  20,000  ft  duration  density  steps  of  +75  and  +100% 
magnitude  with  lower  edges  between  100,000  ft  and  120,000  ft  are  sufficient  to  cause  a 
catastrophic  failure  requiring  more  than  1000  ft/sec  of  propulsive  maneuvering  to  circular¬ 
ize  in  the  desired  orbit.  These  failures  are  caused  because  the  vehicle  enters  the  high  den- 
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sity  region  near  the  bottom  of  the  trajectory.  When  the  onboard  accelerometers  measure 
the  rapid  deceleration  caused  by  the  high  density  pocket  and  feed  this  to  the  control  sys¬ 
tem  the  control  system  responds  by  applying  lift  up  to  decrease  the  energy  rate.  As  the  ve¬ 
hicle  exits  the  high  density  pocket  the  control  system  responds  by  commanding  a  lift  down 
configuration.  But  the  vehicle’s  response  time  is  too  slow,  requiring  thirteen  seconds  to 
perform  a  1 80°  rest  to  rest  roll  maneuver.  By  the  time  the  maneuver  is  complete  the  vehi¬ 
cle  has  moved  higher  in  the  atmosphere  and  no  longer  is  able  to  control  the  trajectory  us¬ 
ing  aerodynamic  forces.  The  vehicle  exits  the  atmosphere  without  properly  depleting  the 
kinetic  energy.  This  behavior  could  also  be  called  a  skipout.  The  same  phenomena  is  ob¬ 
served  for  the  high  magnitude  density  pulses  in  the  nominal  and  high  density  atmospheres, 
though  the  effect  is  less  disastrous. 

The  locations  of  the  density  pulses  which  cause  the  problems  are  higher  in  the  nomi¬ 
nal  atmosphere  than  in  the  low  density  atmosphere,  and  even  higher  still  in  the  high  densi¬ 
ty  atmosphere  than  in  the  nominal  atmosphere.  These  higher  locations  are  because  the 
vehicle’s  initial  configuration  is  lift  up.  The  higher  density  atmospheres  exert  more  aero¬ 
dynamic  force  at  higher  altitudes,  tending  to  decrease  the  vehicle’s  negative  altitude  rate 
earlier  and  increase  the  altitude  at  which  the  vehicle  bottoms  out.  A  density  pulse  which 
perturbs  the  trajectory  near  its  minimum  altitude  must  be  located  higher  in  a  high  density 
atmosphere  than  in  a  low  density  atmosphere. 

One  additional  drawback  to  the  Energy  Controller  is  higher  trajectory  loads  than  the 
algorithms  which  use  the  equilibrium  glide  phase.  The  equilibrium  glide  phase  holds  the 
lift  up  configuration  until  the  trajectory  bottoms  out  in  almost  all  cases.  The  Energy  Con¬ 
troller  will  roll  the  vehicle  from  lift  up  before  the  vehicle  bottoms  out,  allowing  the  vehicle 
to  sink  to  a  lower  minimum  altitude,  producing  higher  peak  aerodynamic  heating  loads, 
higher  maximum  dynamic  pressures,  and  higher  maximum  acceleration  loads.  These 
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higher  trajectory  loads  may  reduce  somewhat  the  advantages  of  aerobraking,  especially  if 
they  require  the  vehicle  to  be  built  heavier  to  withstand  the  acceleration  forces,  or  if  they 
require  additional  heat  shields  or  ablative  materials.  Though  this  study  is  concerned  with 
control  system  robustness,  the  effect  of  a  control  system  on  trajectory  loads  must  be  con¬ 
sidered. 

The  Energy  Controller  shows  some  shortcomings,  especially  with  respect  to  short 
period  density  variations  and  trajectory  loads  which  make  it  unsatisfactory  for  controlling 
a  vehicle  aerobraking  in  the  Martian  atmosphere. 

Mars  Hybrid  Predictor  Corrector  Performance  Results 

The  Mars  Hybrid  Predictor  Corrector  (MHPC)  was  one  of  the  two  best  performing 
algorithms  tested  for  this  series  of  perturbations.  As  discussed  in  Chapter  11  this  control 
algorithm  employs  a  variable  transition  velocity  for  the  switch  from  the  equilibrium  glide 
phase  to  the  predictor  corrector  exit  phase.  Equally  important  is  the  density  estimation 
technique  which  measures  and  records  density  at  discrete  altitude  locations  during  the  en¬ 
try  into  the  atmo.sphere.  Density  during  the  exit  from  the  atmosphere  is  measured  and 
compared  against  that  predicted  using  the  stored  entry  data.  The  result  is  filtered  to  re¬ 
move  high  frequency  noise  and  used  to  bias  the  density  estimate  developed  during  entry. 
The  biased  estimate  is  then  used  to  predict  velocity  loss  for  the  remainder  of  the  trajectory. 
As  may  be  suspected,  this  method  is  extremely  effective  whenever  the  density  profiles  for 
the  inbound  and  outbound  legs  of  the  trajectory  are  the  same. 

The  performance  of  this  controller  is  presented  in  Fig.  28,  Fig.  29,  Fig.  30  and  Fig. 
31.  Again  the  first  three  figures  summarize  the  performance  when  the  density  function  is 
perturbed  with  10,000  and  20,000  ft  duration  square  wave  density  steps.  Again,  the  first. 
Fig.  28,  presents  these  results  when  perturbing  waves  are  injected  into  a  nominal  atmo- 
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Fig.  31  Mars  Hybrid  Predictor  Corrector  Sensitivity  to  Variations  in  Lift  to  Drag 


Ratio  and  Entry  Flight  Path  Angle,  a)  Nominal  Atmosphere;  b)  Low  Density 


Atmosphere;  c)  High  Density  Atmosphere 
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sphere.  Fig.  29  shows  the  results  in  a  low  density  atmosphere  and  Fig.  30  in  a  high  density 
atmosphere.  In  all  of  these  cases  the  controller  was  able  to  guide  the  vehicle  to  the  target¬ 
ed  exit  state  almost  perfectly.  With  this  controller  the  required  AV  following  a  maneuver 
in  a  high  density  atmosphere  is  slightly  less  than  that  in  a  nominal  atmosphere,  which  is 
again  slightly  less  than  that  in  a  low  density  atmosphere.  These  results  are  in  agreement 
with  those  found  using  the  Conjugate  Gradient  optimization  technique  of  Appendix  A. 

The  performance  when  the  vehicle  lift  and  drag  characteristics  are  varied,  and  when 
the  entry  flight  path  angle  are  varied  (Fig.  31)  are  equally  promising.  The  reader  is  cau¬ 
tioned  that  the  results  presented  here  were  all  generated  with  density  functions  which  are 
simply  functions  of  altitude.  The  density  function  for  the  outbound  leg  of  the  trajectory  is 
identical  to  the  density  on  the  inbound  leg.  The  density  estimator  in  this  control  algorithm 
gives  excellent  results  when  the  outbound  density  function  matches  that  measured  while 
inbound,  and  the  control  algorithm  is  able  to  guide  the  vehicle  to  near  perfect  exit  state 
whenever  it  is  supplied  with  a  good  density  function  estimate.  Later,  in  Chapter  5  the  per¬ 
formance  will  be  evaluated  whenever  the  inbound  and  outbound  density  functions  differ. 

Mars  Predictor  Corrector  Performance  Results 

The  Mars  Predictor  Corrector  (MPC)  Control  Algorithm  differs  from  the  MHPC  of 
the  previous  section  only  in  the  density  estimation  technique  employed.  The  MPC  mea¬ 
sures  and  stores  density  every  second  during  the  descent  into  the  atmosphere.  These  den¬ 
sity  measurements  are  then  normalized  using  a  two  stage  exponential  function.  The 
resulting  normalized  data  are  fit  with  a  sixth  order  polynomial  in  altitude.  This  polynomi¬ 
al  is  continually  updated  throughout  the  trajectory  after  each  density  measurement  is  tak¬ 
en.  Again,  the  resulting  density  estimate  is  used  to  compute  the  velocity  loss  yet  to  occur 
due  to  aerodynamic  drag. 
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The  data  generated  during  the  testing  phase  for  the  MPC  were  not  as  good  as  those 
from  the  MHPC.  This  density  estimation  technique  might  be  expected  to  perform  slightly 
worse,  certainly  no  better,  than  the  hybrid  density  estimation  technique  whenever  the  in¬ 
bound  and  outbound  density  functions  are  the  same.  The  strength  of  this  density  estima¬ 
tion  technique  is  expected  to  be  those  cases  when  the  inbound  and  outbound  density 
functions  are  different  (again,  to  be  investigated  in  Chapter  V). 

The  performance  of  this  algorithm  is  presented  in  Fig.  32,  Fig.  33,  Fig.  34  and  Fig. 
35.  The  first  three  figures  present  the  results  of  square  wave  density  pulses,  while  the  last 
figure  shows  the  results  of  varying  the  lift  and  drag  coefficients  and  the  entry  flight  path 
angle. 

The  performance  of  this  algorithm  shown  in  Fig.  32,  Fig.  33,  Fig.  34  and  Fig.  35,  though 
not  as  good  as  that  of  the  MHPC,  is  still  acceptable.  The  worst  performance  noted  here, 
caused  by  a  20,000  ft  duration  +75%  density  pulse  in  the  high  density  atmosphere  located 
between  180,(XX)  and  2(X),000  ft,  required  457  ft/sec  to  attain  the  desired  orbit.  This  algo¬ 
rithm,  when  faced  with  variations  in  entry  flight  path  angle  and  L/D,  produces  practically 
flat  performance  maps  that  are  very  near  the  idealized  optimal  values  calculated  using  the 
method  of  Appendix  A. 


l^V/sec) 


Fig.  35  Mars  Predictor  Corrector  Sensitivity  to  Variations  in  Lift  to  Drag  Ratio 
and  Entry  Flight  Path  Angle,  a)  Nominal  Atmosphere;  b)  Low  Density 
Atmosphere;  c)  High  Density  Atmosphere 
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Lyapunov  Hybrid  Tracking  Controller  Performance  Results 

The  Lyapunov  Hybrid  Tracking  Controller  (LHPC),  derived  in  Chapter  111,  tied  with 
the  MHPC  as  the  two  best  control  algorithms  for  this  test  sequence.  This  controller  em¬ 
ploys  the  equilibrium  glide  phase  for  the  first  part  of  the  trajectory  and  a  Lyapunov  Track¬ 
ing  exit  phase  using  Lyapunov  steepest  descent  techniques  to  steer  the  trajectory  onto  a 
target  path.  The  target  path  selected  is  a  450  ft/sec  constant  altitude  rate  path  which  will 
lead  the  vehicle  to  the  desired  apocenter  altitude.The  transition  velocity  for  switching 
from  the  equilibrium  glide  phase  to  the  LHTC  exit  phase  is  varied  using  the  technique  of 
Chapter  II.  The  hybrid  density  estimation  technique  presented  in  Chapter  II  is  used  to  de¬ 
fine  the  desired  path  and  to  select  the  appropriate  transition  velocity. 

Testing  this  controller  against  the  same  perturbations  considered  earlier  in  this  chap¬ 
ter  produced  excellent  results.  The  results  of  injecting  square  wave  density  pulses  into  the 
nominal  atmosphere,  low  density  atmosphere  and  high  density  atmosphere  are  summa¬ 
rized  in  Fig.  36,  Fig.  37  and  Fig.  38  respectively.  The  results  of  varying  L/D  and  entry 
flight  path  angle  are  presented  in  Fig.  39. 

This  controller  showed  outstanding  results  to  this  test  program  with  practically  no 
sensitivity  to  any  of  these  perturbations.  Again,  however,  the  same  caution  presented  in 
the  MHPC  performance  results  section  should  be  repeated  here:  this  simulated  density  is 
simply  a  function  of  altitude.  The  density  function  for  the  outbound  leg  of  the  trajectory  is 
identical  to  the  density  on  the  inbound  leg.  The  density  e.stimation  technique  employed  in 
this  controller  should  excel  under  this  condition.  The  robustness  with  respect  to  horizontal 
density  variations  must  be  evaluated  to  fairly  generalize  the  evaluation  of  this  (or  any  oth¬ 
er)  guidance  scheme. 


Fig.  39  Lyapunov  Hybrid  Tracking  Controller  Sensitivity  to  Variations  in  Lift  to 
Drag  Ratio  and  Entry  Flight  Path  Angle,  a)  Nominal  Atmosphere;  b)  Low 
Density  Atmosphere;  c)  High  Density  Atmosphere 
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Lyapunov  Tracking  Controller  Performance  Results 

The  LTC  differs  from  the  LHTC  algorithm  only  in  the  density  estimation  technique. 
The  LHTC  employs  the  polynomial  density  estimator  described  in  Chapter  II,  while  the 
LHTC  uses  the  hybrid  density  estimation  technique.  The  performance  of  this  controller 
appears  to  be  slightly  degraded  from  the  performance  of  the  LHTC  at  about  the  same  level 
that  the  performance  of  the  MPC  was  worse  than  that  of  the  MHPC.  The  performance  is 
still  acceptable  and,  as  was  stated  in  the  analysis  of  the  MFC’s  performance,  the  strength 
of  this  density  estimation  technique  is  expected  to  surface  when  the  outbound  and  inbound 
density  functions  differ. 

Fig.  40,  Fig.  41  and  Fig.  42  illustrate  the  results  of  the  square  wave  density  pul«es 
which  perturb  the  nominal,  low  and  high  density  atmosphere,  respectively.  Varying  L/D 
and  entry  flight  path  angle  is  depicted  in  Fig.  43.  The  worst  performance  noted  during 
these  simulations  using  the  LTC  required  464  ft/sec  to  attain  the  desired  orbit.  This  peak 
was  caused  by  a  +75%  density  pulse  perturbing  the  high  density  atmosphere  between 
180,000  and  200,000  ft.  But  again,  even  this  worst  case  is  considered  to  be  acceptable. 

Selection  of  Controllers  to  Proceed 

The  next  stage  of  simulation  was  very  intensive,  requiring  approximately  sixty  hours 
of  computer  time  to  fully  test  each  controller.  In  an  effort  to  limit  this  test  matrix,  only 
those  controllers  which  showed  promise  of  being  able  to  handle  the  perturbations  used  in 
this  chapter  were  to  proceed  to  the  next  phase.  The  original  plan  was  to  select  the  two 
most  promising  controllers,  and  validate  them.  But,  after  analyzing  the  data  presented  in 
this  chapter  four  controllers  were  selected  for  the  next  pha.se  of  testing.  The  four  selected 
were  the  MHPC,  the  MPC,  the  LHTC  and  the  LTC.  The  four  selected  are  actually  two 


Fig.  43  Lyapunov  IVacking  Controller  Sensitivity  to  Variations  in  Lift  to  Drag 
Ratio  and  Entry  Flight  Path  Angle,  a)  Nominal  Atmosphere;  b)  Low  Density 
Atmosphere;  c)  High  Density  Atmosphere 
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control  techniques,  each  employing  two  different  methods  of  density  estimation.  All  four 
of  these  controllers  were  able  to  handle  the  full  range  of  testing  performed  during  this 
phase  without  requiring  more  than  500  ft/sec  to  attain  the  desired  orbit  for  any  perturba¬ 
tion.  The  limitation  of  this  test  sequence  was  that  the  inbound  and  outbound  density  func¬ 
tions  were  always  the  same.  In  the  next  chapter  the  performance  of  these  four  control 
algorithms  will  be  determined  when  the  inbound  and  outbound  density  functions  differ. 
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CHAPTER  V 

DETERMINATION  OF  ROBUSTNESS  LIMITS 

Determination  of  the  robustness  limits  of  the  MHPC,  MPC,  LHTC  and  LTC,  and  se¬ 
lection  of  the  most  robust  algorithm  of  these  four  is  the  goal  of  this  chapter.  Chapter  IV 
showed  that  these  four  algorithms  are  all  capable  of  handling  extreme  variations  in  the  ve¬ 
hicle  L/D,  entry  flight  path  angle  and  in  the  density  function,  provided  the  density  function 
is  a  simple  function  of  altitude.  This  chapter  will  examine  the  effect  of  more  realistic  den¬ 
sity  functions  which  differ  for  the  inbound  and  outbound  legs  of  the  trajectory.  This  will 
be  accomplished  by  again  injecting  square  wave  density  pulses  into  the  density  function, 
but  this  time  the  pulses  will  only  perturb  the  outbound  leg  of  the  trajectory.  In  addition  si¬ 
nusoidal  variations  in  altitude  and  in  vehicle  range  will  be  used  to  perturb  the  density  func¬ 
tion.  The  control  algorithms  will  also  be  tested  using  the  actual  density  profiles  measured 
by  the  Viking  1  and  Viking  2  landers. 

To  determine  the  robustness  limits  success  and  failure  must  first  be  defined.  Because 
the  vehicle  has  not  been  designed  yet,  the  fuel  budget  for  maneuvering  the  vehicle  has  not 
been  defined.  The  definition  of  success  and  failure  used  here  is  somewhat  arbitrary, 
though  it  is  believed  to  be  close  to  the  actual  definition.  Success  is  defined  as  any  aero- 
braking  trajectory  which  requires  500  ft/sec  or  less  of  propulsive  maneuvering  (AV)  to  at¬ 
tain  the  desired  final  orbit.  As  in  Chapter  IV  AV  is  computed  with  three  components,  one 
propulsive  maneuver  applied  at  the  atmospheric  interface  in  the  direction  of  the  velocity 
vector  to  correct  the  apocenter  altitude,  a  second  at  apocenter  to  raise  pericenter  and  a 
third  to  correct  any  plane  error.  Because  of  the  lack  of  a  firm  definition  of  vehicle  charac¬ 
teristics  a  grey  area  has  been  defined.  The  grey  area  is  any  trajectory  which  requires  be¬ 
tween  500  and  1 ,000  ft/.sec.  Any  trajectory  which  requires  between  500  and  1000  ft/sec  to 
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attain  the  desired  orbit  will  be  referred  to  as  a  soft  failure.  Any  vehicle  design  should  cer¬ 
tainly  carry  enough  fuel  to  perform  500  ft/sec  of  maneuvering  to  attain  the  desired  orbit. 
Carrying  additional  fuel,  however,  to  correct  for  a  soft  failure  of  the  aerobraking  control 
system  reduces  the  advantage  of  aerobraking.  The  trajectories  which  terminate  with  soft 
failures  would  still  be  able  to  attain  orbit,  just  not  the  desired  orbit,  using  500  ft/sec  of  pro¬ 
pulsive  maneuvering.  This  anomaly  may  result  in  some  mission  degradation,  but  not  a 
complete  mission  failure.  A  hard  failure  is  defined  to  be  any  trajectory  which  requires 
1 ,000  ft/sec  or  more  of  AV  to  attain  the  desired  orbit.  It  includes  any  trajectories  which  fail 
to  exit  the  atmosphere.  By  this  definition,  all  four  controllers  considered  in  this  chapter 
succeeded  in  all  of  the  simulations  performed  in  Chapter  IV. 

Outbound  Leg  Square  Wave  Density  Pulses 

The  robustness  test  procedure  begins  by  using  square  wave  density  pulses  which  per¬ 
turb  the  density  function  of  the  outbound  leg  only.  The  density  during  the  descent  into  the 
atmosphere  is  either  the  nominal  or  a  MARS-GRAM  generated  low  or  high  density  atmo¬ 
sphere  model.  After  the  altitude  rate  becomes  positive,  a  square  wave  density  pulse,  simi¬ 
lar  to  those  employed  in  Chapter  IV  is  used  to  perturb  either  a  10,000  or  a  20,000  ft 
altitude  band  of  the  atmosphere.  These  pulses  multiply  the  density  predicted  by  the  atmo¬ 
sphere  model  by  0.5,  0.75,  1.25,  1.5,  1.75  or  2.0  within  the  perturbed  altitude  band.  The 
pulses  are  again  referred  to  as  -50%,  -25%,  +25%,  +50%,  +75%  and  +100%  magnitude 
density  pulses  respectively.  As  in  Chapter  IV  the  pulses  are  moved  in  10,000  ft  altitude  in¬ 
tervals,  with  the  lower  edge  of  the  density  pulse  located  between  100,000  and  290,000  ft. 
The  performance  is  presented  in  Fig.  44  through  Fig.  55  with  AV  plotted  along  the  vertical 
axis  while  the  magnitude  of  the  pulse  and  the  location  of  the  lower  edge  are  plotted  on  the 
other  two  axis. 
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MHPC  Performance 

Performance  of  the  MHPC  in  the  nominal,  low  and  high  density  atmospheres  when 
the  density  of  the  outbound  leg  of  the  trajectory  is  perturbed  by  square  wave  density  pulses 
is  presented  in  Fig.  44,  Fig.  45  and  Fig.  46  respectively.  In  the  first  plot  of  each  figure  the 
pulse  perturbs  a  10,000  ft  altitude  band,  while  in  the  second  plot  the  pulse  perturbs  a 
20,000  ft  altitude  band. 

The  MHPC  produced  many  soft  failures  during  this  test  sequence  but  no  hard  fail¬ 
ures  were  recorded.  The  20,000  ft  duration  pulses  produce  worse  performance  than  the 
10,000  ft  duration  pulses  in  almost  all  cases.  The  MHPC  is  very  sensitive  to  large  magni¬ 
tude  (+50%  and  +75%)  density  pulses  located  below  180,0(X)  ft  in  the  nominal  atmo¬ 
sphere,  and  150,0(X)  or  200,000  ft  in  the  low  or  high  density  atmosphere  respectively. 
There  is  also  a  region  of  sensitivity  caused  by  the  -50%  20,000  ft  density  pulses.  These, 
though,  are  located  at  slightly  higher  altitudes  and  are  not  as  severe  as  those  caused  by  the 
large  magnitude  pulses.  In  all  of  these  plots  there  is  a  region  at  extremely  low  altitudes, 
where  the  pulses  have  minimal  or  no  effect.  This  robust  region  occurs  because  these  puls¬ 
es  are  either  located  below  the  minimum  altitude  of  the  trajectory  and  the  satellite  never 
flies  in  the  perturbed  atmosphere,  or  they  are  very  near  the  minimum  altitude  of  the  trajec¬ 
tory  and  the  satellite  does  not  spend  much  time  in  the  perturtted  atmosphere. 

There  are  two  primary  failure  modes  for  these  trajectories.  When  the  large  magni¬ 
tude  density  pulse  perturbs  the  atmosphere  in  the  altitude  region  where  the  satellite  is  in 
the  equilibrium  glide  phase,  the  density  filter  is  fooled  into  believing  the  entire  atmosphere 
has  higher  density  than  that  measured  during  the  descent.  The  effect  is  to  initiate  the  exit 
phase  early,  and  predict  a  relatively  high  altitude  rate  for  the  exit  phase.  When  the  vehicle 
moves  out  of  the  high  density  region  there  is  a  time  lag  before  the  density  filter  records  the 
change.  By  the  time  the  controller  responds  the  vehicle  has  moved  even  higher,  and  the 
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vehicle  has  less  control  authority.  The  result  is  the  vehicle  leaves  the  atmosphere  with  too 
much  energy,  and  overshoots  the  desired  apocenter  altitude.  The  second  failure  mode  is 
caused  by  higher  altitude  density  pulses.  All  of  the  failures  caused  by  the  -50%  pulses  ex¬ 
hibited  this  failure  mode.  The  vehicle  flies  the  equilibrium  glide  phase  in  the  unperturbed 
atmosphere.  After  the  vehicle  initiates  the  exit  phase  it  encounters  the  perturbed  atmo¬ 
sphere.  The  large  magnitude  density  pulses  dissipate  more  energy  than  predicted  resulting 
in  a  steeper  exit  phase  than  desired,  and  in  some  cases  a  lower  apocenter  than  desired. 
Conversely,  the  small  magnitude  density  pulses  cause  the  vehicle  to  lose  less  energy  than 
predicted  and  result  in  an  apocenter  altitude  higher  than  desired. 

MPC  Performance 

The  MPC  definitely  has  better  performance  than  the  MHPC  under  these  conditions. 
Again,  the  performance  is  presented  in  three  figures.  Fig.  47,  Fig.  48  and  Fig.  49  with  the 
first  figure  showing  results  from  the  nominal  atmosphere,  the  second  from  the  low  density 
atmosphere  and  the  third  from  the  high  density  atmosphere. 

The  10,000  ft  density  pulses  had  almost  no  effect  on  the  performance  of  this  control 
algorithm.  Even  the  20,000  ft  pulses  produced  reasonably  good  results.  There  was  only 
one  soft  failure  noted  during  this  lest  sequence  and  two  very  near  failures  for  the  MPC. 
All  three  of  these  events  were  caused  by  20,000  ft  +100%  density  pulses  perturbing  the 
low  density  atmosphere.  The  pulse  between  120,000  and  140,000  ft  required  a  AV  of  584 
ft/sec,  while  the  pulse  10,000  ft  higher  (between  130,000  and  150,(XX)  f!)  required  498ft/ 
sec.  The  pulse  between  140,(XX)  and  160,000  ft  required  just  458  ft/sec,  but  the  one  be¬ 
tween  150,(XX)  and  170,(XX)  ft  required  492  fi/sec.  These  high  AVs  were  caused  by  the 
same  failure  modes  as  described  above  with  the  pulses  with  lower  edges  at  120,(X)0  and 
130,0(X)  ft  causing  the  density  estimator  to  overreact  and  force  the  vehicle  to  exit  with  too 
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much  energy,  while  the  pulse  at  150,000  ft  causes  the  vehicle  to  lose  more  energy  than 
planned.  Overall,  though,  the  polynomial  density  estimator  used  in  the  MPC  showed  im¬ 
provement  over  the  h  ’brid  density  estimator  of  the  MHPC.  As  suggested  in  chapter  IV, 
the  promise  of  this  density  estimation  technique  is  realized  when  the  inbound  and  out¬ 
bound  density  functions  are  different. 

LHTC  Performance 

Performance  of  the  LHTC  was  mixed.  For  the  majority  of  density  perturbations  cal¬ 
culated  here  this  controller  performed  better  than  either  the  MHPC  or  the  MPC.  Yet  there 
were  a  few  isolated  instances  where  the  controller  performed  extremely  poorly.  The  con¬ 
troller  even  produced  two  hard  failures.  Presentation  of  this  controller’s  |)erformance  fol¬ 
lows  the  same  format  as  before  with  the  nominal  atmosphere  results  in  Fig.  50,  the  low 
density  atmosphere  results  in  Fig.  51  and  the  high  density  atmosphere  results  in  Fig.  52. 

One  noteworthy  aspect  of  the  Lyapunov  Tracking  exit  phase  is  that  it  almost  always 
commands  either  full  lift  up,  or  full  lift  down.  When  the  vehicle  is  right  on  the  desired 
path  the  commanded  bank  angle  will  chatter  between  ±15°  and  ±165°  (commanded  bank 
angles  less  than  1 5°  or  greater  than  165°  are  allowed  only  when  the  orbit  plane  error  is  less 
than  .03°).  Of  course,  the  vehicle  roll  rate  and  roll  acceleration  limits  prevent  the  vehicle 
from  oscillating  too  wildly.  But,  when  the  vehicle  is  not  on  the  desired  path  the  control 
system  will  command  near  full  lift  up,  or  full  lift  down  to  approach  the  trajectory.  This 
feature  allows  the  vehicle  to  respond  more  quickly  than  it  does  for  the  predictor  corrector 
algorithms  to  pull  the  vehicle  back  onto  the  desired  path.  However,  when  the  desiied  path 
is  computed  poorly  because  of  a  poor  density  estimate,  the  control  system  still  responds  by 
commanding  full  lift  to  approach  the  computed  path  as  rapidly  as  possible.  This  controller 
suffers  from  the  same  problems  with  the  density  estimator  as  the  MHPC.  This  phenomena 
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caused  the  two  hard  failures  seen  in  Fig.  50b  and  Fig.  52b.  The  onboard  accelerometer 
measurements  had  been  recording  high  drag  measurements  while  the  vehicle  was  deceler¬ 
ating  in  the  region  of  high  density  caused  by  the  density  pulse.  The  density  filter  had  pre¬ 
dicted  that,  because  of  the  high  drag  inferred  density  measurements,  the  remainder  of  the 
atmosphere  would  also  be  higher  density  than  that  recorded  during  the  descent.  The  effect 
is  to  predict  higher  energy  loss  due  to  drag  than  will  actually  occur.  This  prediction  causes 
the  control  system  to  initiate  the  exit  phase  earlier  than  desired  and  to  plot  a  path  which 
climbs  out  of  the  atmosphere  at  relatively  high  speed.  The  Lyapunov  optimal  control  solu¬ 
tion  is  to  pull  onto  this  path  as  rapidly  as  possible.  By  the  time  the  satellite  has  moved  out 
of  the  high  density  region,  and  the  density  filter  has  recovered,  the  satellite  has  moved  too 
high,  at  a  velocity  which  is  too  high  too  allow  recovery.  The  result  was  a  post  aerobraking 
apocenter  altitude  of  836  nm  for  the  hard  failure  in  Fig.  50b  and  1,165  nm  in  Fig.  52b. 
Again,  the  target  apocenter  is  270  nm. 

The  Lyapunov  Tracking  exit  phase,  however,  seems  be  able  to  cope  with  these  densi¬ 
ty  estimation  problems  better,  in  most  instances,  than  the  predictor  corrector  algorithms. 
The  rapid  response  of  the  vehicle,  due  to  the  nature  of  the  Lyapunov  control  system, 
though  it  caused  the  two  hard  failures  discussed  above,  was  usually  advantageous.  As  ac¬ 
celerometer  measurements  are  taken,  and  the  density  filter  is  continually  updated,  the  de¬ 
sired  path  varies.  The  rapid  response  of  the  Lyapunov  control  law  helps  track  this  moving 
path  as  long  as  the  vehicle  has  enough  aerodynamic  control  authority  to  respond.  In  the 
LTC  controller  results  which  follow,  the  effect  of  combining  the  polynomial  density  esti¬ 
mation  technique  with  the  fast  response  of  the  Lyapunov  control  scheme  is  investigated. 


8 


Densil 


Ill 


LTC  Performance 

The  LTC  performance  is  presented  in  Fig.  53,  Fig.  54  and  Fig.  55.  The  10,000  ft 
density  pulses  had  almost  negligible  effect  on  the  performance  of  this  control  algorithm  as 
was  the  case  with  the  MFC.  Also,  the  20,000  ft  pulses  produced  reasonably  good  results. 
There  were  only  three  soft  failures  noted  during  this  test  sequence  for  the  LTC.  All  three 
of  the  failures  were  caused  by  the  same  20,000  ft  +100%  density  pulses  perturbing  the  low 
density  atmosphere  which  caused  the  soft  failure  and  the  two  other  near  failures  in  the 
MPC  performance.  For  the  LTC  the  pulse  between  120,000  and  140,000  ft  required  a  AV 
of  542  ft/sec,  while  the  pulse  10,000  ft  higher  (between  130,000  and  150,(XX)  ft)  required 
505  ft/sec.  The  pulse  between  140,(X)0  and  160,000  ft  did  not  result  in  a  failure  requiring 
493  ft/sec,  but  the  one  between  150, (XX)  and  170,000  ft  did,  requiring  533  ft/sec.  These 
failures  were  caused  by  the  same  failure  modes  as  described  earlier  with  the  pulses  locat¬ 
ed  at  120,000  and  130,000  ft  causing  the  density  estimator  to  overreact  forcing  the  vehicle 
to  exit  with  too  much  energy,  while  the  pulse  at  150,0(X)  ft  causes  the  vehicle  to  lose  more 
energy  than  planned. 

Overall,  though,  the  polynomial  density  estimator  u.sed  in  combination  with  the  Ly¬ 
apunov  control  scheme  shows  excellent  performance.  The  performance  of  this  control  al¬ 
gorithm  during  this  test  sequence  very  nearly  paralleled  that  of  the  MPC. 

Sinusoidal  Density  Variations 

The  next  testing  sequence  involves  perturbing  the  density  function  with  sine  waves. 
Sine  waves  in  altitude  and  sine  waves  in  range  were  used.  These  sine  waves  were  varied  in 
amplitude  wavelength  (X.)  and  phase  angle  (({»).  The  sine  wave  perturbations  in  alti-  , 
tude  took  the  form 


l^\,/sec) 


115 
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while  those  in  range  took  the  form 
P  ~  Pmodel 
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The  range  of  amplitudes  used  included  0.1,  0.25  and  0.5  for  both  forms  of  perturba¬ 
tions.  The  wavelengths  selected  for  the  altitude  variations  included  1000,  2000,  5000, 
10,000,  20,000,  50,000,  100,000,  200,000  and  500,000  feet.  For  the  variations  in  range 
the  wavelengths  selected  included  1,5,  10,  20,  50,  100,  200,  500,  1000,  2000  and  5000 
nautical  miles.  In  both  cases,  the  phase  angles  included  zero  through  ^  in  ^  increments. 


The  sinusoidal  variations  with  amplitude  of  0.1  appears  to  be  very  much  in  line  with 
the  actual  density  profiles  measured  by  Viking  1  and  Viking  2  landers  during  their  descent 
through  the  Martian  atmosphere  (Fig.  7  and  Fig.  8).  All  of  the  results  with  an  amplitude  of 
0. 1  for  both  forms  of  the  sinusoidal  variations  and  all  four  controllers  examined  in  this 
chapter  were  successful.  The  highest  AV  required  490  ft/sec,  but  the  vast  majority  of  the 
trajectories  (over  99%)  required  less  than  400  ft/sec.  Only  13  of  the  1920  trajectories  test¬ 
ed  with  a  0.1  amplitude  sine  wave  density  variation  required  more  than  400  ft/sec  of  AV. 
Likewise,  the  results  generated  using  25%  and  50%  amplitude  sine  waves  in  altitude  are 
almost  as  benign  as  the  10%  results.  All  of  the  trajectories  which  used  25%  amplitude 
sine  waves  in  altitude  were  successful.  Of  the  864  trajectories  checked  using  50%  ampli¬ 
tude  sine  waves  in  altitude  none  resulted  in  hard  failures  and  only  8  produced  soft  failures. 
Of  these,  only  three  required  more  than  600  ft/sec  with  the  worst  requiring  732  ft/sec.  A 
complete  breakdown  of  these  failures  is  presented  in  Table  2.  Because  these  results  gener¬ 
ated  using  10%  sine  waves  in  altitude  and  range  and  25%  sine  waves  in  altitude  were  all 
successful,  and  the  eight  soft  failures  generated  using  50%  sine  waves  in  altitude  are  ade¬ 
quately  described  in  Table  2  they  will  not  be  presented  graphically.  It  is  interesting  to  note 
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that  three  of  the  four  soft  failures  which  occurred  with  controllers  using  the  hybrid  density 
estimation  technique  had  wavelengths  of  1000  ft,  while  the  fourth  had  a  wavelength  of 
5000  ft.  The  1000  ft  wavelength  sine  waves  in  altitude  seem  to  be  corrupting  the  stored 
density  data  used  in  the  density  estimation  process.  This  data  is  stored  at  1000  ft  altitude 
intervals.  Though  the  density  filter  should  be  able  to  compensate  for  this,  it  does  not  ap¬ 
pear  to  do  so  well  enough  to  prevent  these  failures.  Three  of  the  four  failures  which  oc¬ 
curred  with  controllers  employing  the  polynomial  density  estimation  technique  had 
wavelengths  of  20,000  and  50,000  ft.  Shorter  wavelengths  tend  to  have  a  cancelling  ef¬ 
fect,  with  the  additional  drag  of  high  density  regions  being  offset  by  the  lower  drag  of  low 
density  regions.  Longer  wavelengths  are  easy  for  the  sixth  order  polynomial  to  follow, 
provided  there  are  no  more  than  five  extremes  in  the  density  function.  The  problem  with 
the  20,000  and  50,000  ft  wavelength  sine  waves  is  they  do  not  oscillate  fast  enough  to  can¬ 
cel  high  density  regions  against  low  density  regions,  yet  they  still  have  six  to  fifteen  com¬ 
plete  sine  waves,  with  twelve  to  thirty  density  extremes  in  the  aerobraking  region;  more 
than  a  sixth  order  polynomial  can  follow.  The  final  failure  was  caused  by  an  excessive  or¬ 
bit  plane  error  (wedge  angle)  at  exit. 

25%  and  50%  Sine  Waves  in  Range 

The  25%  amplitude  sine  wave  density  perturbations,  which  use  vehicle  range  from 
entry  as  the  argument  to  the  sine  function,  are  probably  the  truest  measure  used  here  to  test 
controller  robustness  in  the  presence  of  a  realistic  worst  case  Martian  atmosphere.  These 
perturbations  are  of  somewhat  higher  magnitude  than  the  perturbations  measured  by  the 
Viking  1  and  Viking  2  landers,  but,  most  probably,  the  Viking  1  and  Viking  2  landers  did 
not  sample  the  worst  case  atmospheric  perturbations.  Though  the  amplitude  of  the  per¬ 
turbing  sine  wave  is  increased  to  50%  for  the  test  sequence,  the  probability  is  extremely 
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low  that  the  Martian  atmosphere  would  ever  experience  high  frequency  oscillations  in 
density  with  this  large  of  amplitude. 

The  performance  of  the  MHPC  when  tested  against  the  25%  amplitude  perturbations 
is  presented  in  Fig.  56a,  Fig.  57a  and  Fig.  58a  for  the  nominal,  low  density  and  high  densi¬ 
ty  atmospheres  respectively.  Fig.  56b,  Fig.  57b,  and  Fig.  58b  present  the  results  when  the 
amplitude  of  the  perturbing  sine  wave  is  increased  to  50%.  Similarly,  Fig.  59,  Fig.  60  and 
Fig.  61  present  the  results  for  the  MFC,  while  the  results  for  the  LHTC  is  presented  in 
Fig.  62,  Fig.  63  and  Fig.  64  and  the  LTC  results  are  presented  in  Fig.  65,  Fig.  66  and  Fig. 
67. 


The  25%  amplitude  perturbations  were  significant  enough  to  cause  problems  for 
some  of  the  trajectories.  Though  they  did  not  induce  any  hard  failures,  there  were  many 
soft  failures.  The  50%  amplitude  perturbations  were  severe  enough  to  cause  several  hard 
failures  for  all  of  the  controllers  except  the  LTC.  The  25%  and  the  50%  amplitude  sine 
waves  were  each  used  to  simulate  264  perturbed  atmospheres  for  each  controller  (11 
wavelengths  x  8  phase  angles  x  3  ba.se  atmospheres).  Of  these  264  trajectories  tested  with 
the  MHPC  with  the  25%  amplitude  variation  six  trajectories  resulted  in  soft  failures.  Six 
trajectories  also  resulted  in  soft  failures  when  the  MPC  controller  was  used,  though  they 
were  not  the  same  six  perturbations.  The  LHTC  had  four  si'ft  failures  while  the  LTC  only 
had  two.  When  the  amplitude  of  the  perturbing  sine  wave  was  increased  to  50%  the 
MHPC  had  fourteen  hard  failures,  the  MPC  had  ten  and  the  LHTC  had  fourteen.  These 
three  controllers  also  experienced  many  soft  failures  during  these  simulations.  The  LTC 
did  not  result  in  any  hard  failures,  but  it  did  produce  twenty  nine  .soft  failures. 

All  of  these  failures  were  the  result  of  exit  phase  failures,  which  are  in  turn  attribut¬ 
able  to  density  estimation  problems.  The  equilibrium  glide  phase  was  robust  enough  to 
l.»ep  the  vehicle  in  the  atmosphere  and  prevent  a  skip  out  for  all  of  these  trajectories. 
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None  of  the  trajectories  failed  to  exit  the  atmosphere,  although  some  of  them  barely  did. 
The  problem  with  all  of  the  failures  centered  around  the  inability  of  the  density  estimation 
technique  to  adequately  predict  the  density  function  and  the  amount  of  urug  thus  expected 
by  the  controller  for  the  duration  of  the  trajectory.  Both  density  estimation  techniques  ap¬ 
propriately  ignored  the  high  frequency  density  variations  (those  wavelength  less  than 
10  nm).  These  oscillations  occur  so  quickly  that  the  high  and  low  density  regions  have  a 
cancelling  effect. 

The  hybrid  density  estimator  shows  increased  sensitivity  to  wavelengths  of  20  to  200 
nm.  The  polynomial  density  estimator,  on  the  other  hand  handles  these  wavelengths  very 
well.  It  is  the  500  to  2000  nm  wavelengths  which  produce  problems  for  this  estimator. 
These  sensitivities  to  different  wavelengths  are  easy  to  understand.  The  hybrid  density  es¬ 
timation  technique  uses  the  density  filter  to  adjust  it’s  estimate  for  the  entire  atmosphere 
based  on  the  current  density  measurements.  The  long  wavelength  sine  waves  have  the 
same  effect  as  a  slowly  increasing  or  decreasing  density  bias  during  the  trajectory.  The 
density  filter  of  the  hybrid  density  estimator  is  able  to  sense  this  slow  drift  and  appropri¬ 
ately  adjust  the  measurements  taken  during  descent  to  comjiensate  for  the  drift.  The  wave¬ 
lengths  which  gives  the  hybrid  density  estimator  trouble  are  those  which  perturb  a  portion 
of  the  atmosphere  and  then  reverse  that  perturbation  fast  enough  to  confuse  the  density  fil¬ 
ter  but  not  fast  enough  to  have  a  cancelling  effect.  The  polynomial  density  estimator,  on 
the  other  hand,  fits  the  sixth  order  polynomial  in  altitude  to  the  normalized  density  func¬ 
tion.  This  density  estimation  technique  remembers  the  density  which  was  measured  at  the 
various  altitude  intervals.  It  takes  the  most  recent  density  measurement  and  adds  this  to 
the  knowledge  base  and  fits  a  smooth  polynomial  curve  through  the  data.  When  the  local 
density  is  biased,  but  then  that  bias  reverses  later  in  the  trajectory,  as  it  does  when  the  in¬ 
termediate  wavelength  sine  waves  perturb  the  atmosphere,  this  density  estimation  tech¬ 
nique  excels.  But,  when  the  density  function  is  monotonically  increasing  or  decreasing 
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during  the  trajectory,  as  is  the  case  for  the  longer  wavelength  sine  waves,  this  estimation 
technique  does  not  respond  fast  enough.  An  attempt  to  place  more  weight  on  the  most  re¬ 
cent  data  would  seem  to  help  this  process,  but  attempts  to  do  so  made  the  oldest  data  obso¬ 
lete;  that  is,  the  higher  altitude  densities,  with  the  density  estimator  sometimes  missing  the 
density  at  exit  by  an  order  of  magnitude  or  more.  Clearly,  this  area  deserves  further  study. 

Overall,  the  Lyapunov  control  scheme  performed  better  than  the  predictor  corrector. 
The  rapid  response  of  the  Lyapunov  tracking  exit  phase  was  able  to  compensate  for  the 
slowly  developing  density  estimates.  The  polynomial  density  estimator  also  performed 
better  than  the  hybrid  density  estimator,  as  can  be  seen  in  the  LTC  results.  The  LTC  kept 
AV  below  500  ft/sec  for  all  but  two  of  the  25%  amplitude  sine  wave  perturbed  atmo¬ 
spheres,  and  those  two  only  required  509  and  577  ft/sec.  Additionally,  the  50%  amplitude 
trajectories  were  all  completed  with  AV  below  1000  ft/sec.  The  LTC  was  also  able  to  cope 
with  the  square  wave  density  pulses,  both  those  presented  in  Chapter  IV  which  perturbed 
the  entire  atmosphere  and  those  of  this  chapter  which  only  effect  the  outbound  leg  of  the 
trajectory.  Since  the  LTC  required  less  than  500  ft/sec  for  all  of  the  trajectories  tested  in 
Chapter  IV,  and  responded  better  than  any  of  the  other  controllers  to  all  the  robustness 
tests  of  this  chapter  the  LTC  is  selected  as  the  most  robust  aerobraking  controller  exam¬ 
ined. 
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CHAPTER  VI 

CONCLUSIONS  AND  RECOMMENDATIONS 


Conclusions 

The  Analytic  Predictor  Corrector  algorithm  selected  as  the  control  algorithm  for  the 
AFE  is  generally  a  robust  control  algorithm,  especially  with  respect  to  large  scale  density 
variations.  The  algorithm  is  fairly  robust  to  short  period  density  variations,  but  does  dem¬ 
onstrate  a  definite  sensitivity  to  variations  in  the  entry  flight  path  angle  and  vehicle  lift  and 
drag  coefficients.  These  sensitivities  are  due,  in  large  part,  to  the  fixed  transition  velocity 
employed  to  switch  the  control  algorithm  from  the  entry  phase  to  the  exit  phase  and  the 

rather  simplistic  density  estimation  scheme  used.  It  is  necessary  to  increase  the  K-  term 

*1 

in  the  equilibrium  glide  phase  to  prevent  rapid  large  scale  density  variations  from  causing 
a  premature  exit  from  the  Martian  atmosphere. 

The  Energy  Controller  was  slightly  more  robust  than  the  APC  to  variations  in  the  en¬ 
try  flight  path  angle,  and  vehicle  lift  and  drag  coefficients.  It  was  also  robust  to  large  scale 
density  variations.  However,  short  period  density  variations  were  murderous  to  this  con¬ 
trol  algorithm  and  the  increased  trajectory  loads  caused  by  the  EC  led  to  its  early  dismissal 
from  the  list  of  potential  control  algorithms. 

The  Numerical  Gradient  technique,  and  then  the  Conjugate  Gradient  technique  were 
used  to  compute  idealized  optimal  (minimum  AVj  trajectories.  It  was  hoped  that  these 
methods  could  be  adapted  as  an  on-board  control  algorithm.  But,  these  algorithms  require 
about  two  orders  of  magnitude  more  computational  time  than  the  APC  or  EC  to  generate  a 
solution.  Additionally,  the  optimization  technique  assumes  all  pertinent  density  and  vehi¬ 
cle  lift  and  drag  characteristics  are  known  preci.sely.  The  trajectories  produced  by  these 
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optimization  techniques  fly  the  exit  phase  using  the  full  lift  available  to  remain  in  the  at¬ 
mosphere.  Any  decrease  in  density  from  that  modeled  in  the  optimization  process  allows 
the  vehicle  to  exit  early  with  too  much  velocity.  These  algorithms,  with  the  current  perfor¬ 
mance  index,  would  not  produce  robust  trajectories  even  if  they  were  able  to  compute  a  so¬ 
lution  fast  enough  to  be  used  in  real  time  to  control  the  satellite.  A  more  general 
performance  index  which  seeks  to  minimize  AV  while  retaining  robustness  and  also  reduc¬ 
ing  control  activity  should  be  sought  if  these  techniques  are  to  become  practical. 

The  modifications  proposed  to  the  APC  to  produce  the  MHPC  and  MPC  convert  that 
algorithm  into  a  robust  control  algorithm  capable  of  guiding  the  aerobraking  trajectory  to 
near  minimum  AV  exit  state  for  most  of  the  perturbations  considered.  As  mentioned  be¬ 
fore,  it  was  necessary  to  increase  K-  for  the  equilibrium  glide  phase  to  prevent  a  prema- 

fl 

ture  exit  from  the  atmosphere.  But  in  addition,  the  change  to  the  more  computationally 
straight  forward  and  efficient  exit  phase  combined  with  the  better  density  estimation  tech¬ 
niques  and  the  variable  transition  velocity  made  significant  headway  in  improving  the  ro¬ 
bustness  of  the  control  algorithms.  Between  the  MHPC  and  the  MPC,  the  MPC  responded 
better  overall  to  the  perturbations  examined  here.  There  were  two  areas  where  the  MHPC 
did  slightly  better  than  the  MPC.  The  first  was  the  situations  when  density  is  simply  a 
function  of  altitude  and  the  entry  and  exit  density  functions  were  identical.  This  situation 
is  probably  rather  unrealistic  and  the  MPC  was  still  able  to  handle  these  situations  well 
(though  not  as  well  as  the  MHPC),  without  producing  any  failures.  The  second  area  was 
when  the  large  amplitude  sinusoidal  variations,  which  used  range  from  entry  as  the  argu¬ 
ment  to  the  sine  function,  and  had  wavelengths  between  500  and  2000  nm.  This  area  is 
still  a  concern  and  leads  to  several  of  the  recommendations  below.  Overall,  however,  the 
MPC  reacted  more  appropriately  to  realistic  perturbations  than  did  the  MHPC. 
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The  Lyapunov  Steepest  Descent  Control  algorithm  was  implemented,  but  its  inabili¬ 
ty  to  compensate  for  varied  energy  depletion  rates  due  to  density  variations,  variations  in 
entry  flight  path  angle  or  vehicle  drag  coefficient  made  this  algorithm  unusable.  However, 
when  the  Lyapunov  control  algorithm  was  recast  as  a  tracking  controller  designed  to  fol¬ 
low  a  reference  trajectory,  it  showed  much  more  promise.  The  algorithm  still  had  trouble 
exiting  with  just  the  right  amount  of  energy  to  target  the  desired  apocenter  altitude  and 
produced  peak  trajectory  loads  higher  than  those  of  the  predictor  corrector  algorithms.  To 
cure  the  first  ailment  a  scheme  to  vary  the  gain  values  in  the  Lyapunov  function  was  devel¬ 
oped,  while  the  second  was  fixed  by  employing  the  equilibrium  glide  entry  phase  and  us¬ 
ing  the  Lyapunov  Tracking  Algorithm  as  an  exit  phase. 

With  the  two  density  estimation  techniques  developed  for  the  MHPC  and  the  MPC 
used  to  define  the  reference  trajectory,  and  the  transition  velocity  from  entry  to  exit  phase 
computed  as  for  the  predictor  correctors,  the  LHTC  and  LTC  performed  extremely  well. 
The  performance  of  the  LHTC  and  LTC  essentially  mirrored  that  of  the  MHPC  and  MPC 
respectively.  Generally,  the  strengths  of  the  MHPC  turned  out  to  be  the  strong  points  of 
the  LHTC,  while  they  shared  common  weaknesses  as  well.  Likewise,  perturbations  which 
caused  problems  for  the  MPC  were  also  likely  to  cause  problems  for  the  LTC.  In  most 
cases,  the  problems  were  initiated  because  the  density  estimation  technique  was  unable  to 
follow  a  specific  perturbation.  The  Lyapunov  tracking  algorithm,  with  it’s  more  rapid  re¬ 
sponse,  was  able  to  compensate  better  and  produce  exit  states  which  required  less  AV  than 
the  predictor  correctors.  There  were  a  few  notable  exceptions  where  the  rapid  response 
moved  the  vehicle  into  a  less  dense  region  too  rapidly  resulting  in  loss  of  control  authority 
and  an  exit  state  with  too  much  energy.  But,  predominantly,  the  Lyapunov  trackers  per¬ 
formed  better  than  the  predictor  correctors.  As  in  the  predictor  corrector  analysis  the  poly¬ 
nomial  density  estimation  technique  worked  better  than  the  hybrid  density  estimation 
technique.  Overall,  the  LTC  performed  belter  than  the  LHTC,  MPC  or  MHPC  and  is  the 
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recommended  control  algorithm  for  performing  an  interpl  inetary  aerobraking  maneuver  at 
Mars. 


Recommendations 

Based  on  these  conclusions,  the  following  recommendations  are  made: 

1)  Robustness  to  density  variations  should  be  a  prime  issue  in  selecting  the  control  al¬ 
gorithm  for  the  aerobraking  phase  of  the  MRSR.  This  characteristic  must  be  consid¬ 
ered  along  with  decisions  such  as  entry  velocity,  vehicle  lift  requirements,  ballistic 
coefficient,  or  navigational  accuracy  requirements. 

2)  The  expected  wavelengths  and  maximum  amplitude  of  the  short  period  density  oscil¬ 
lations  in  the  Martian  atmosphere  should  be  characterized.  The  nature  of  these  short 
period  oscillations  should  be  determined.  It  would  be  beneficial  in  designing  a  den¬ 
sity  estimation  technique  to  know  if  the  short  period  density  wave  structure  is  prima¬ 
rily  horizontal  or  vertical  in  nature,  or  a  predominantly  time  varying  function. 

3)  Once  the  frequency  of  the  expected  density  variations  is  determined,  the  density  esti¬ 
mation  technique  employed  in  the  aerobraking  control  system  should  be  tuned  to  re¬ 
spond  to  the  most  likely  frequencies  which  may  perturb  the  trajectory,  while 
ignoring  those  which  have  minimal  effect  on  the  trajectory. 

4)  A  higher  order  density  estimator,  perhaps  using  Tschebechev  polynomials  or  Leg¬ 
endre  polynomials  to  bypass  the  numerical  difficulties  of  a  higher  order  polynomial 
in  altitude  should  be  examined.  It  may  also  be  desirable  to  fit  a  second  function,  in 
terms  of  arc  length,  or  time,  or  range  to  the  density  function,  especially  if  a  monoton- 
ically  increasing  or  decreasing  density  function  is  predicted. 
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5)  The  density  estimation  technique  should  be  adjusted  to  use  all  available  knowledge 
of  the  Martian  atmosphere,  including  any  knowledge  of  dust  storms,  solar  flares  or 
knowledge  of  the  solar  heating  of  the  atmosphere  along  the  intended  trajectory. 

6)  The  LTC  should  be  tested  using  higher  entry  velocities,  different  vehicle  lift  and  drag 
characteristics  or  ballistic  coefficient  as  well  as  different  target  orbits  to  determine  its 
suitability  for  controlling  some  of  the  other  mission  scenarios  proposed  for  MRSR, 
including  the  fast  trip  manned  precursor  mission.  Also,  trading  off  nominal  perfor¬ 
mance  for  robustness  by  varying  the  exit  phase  altitude  rate  should  be  studied. 

7)  A  statistical  method  of  evaluating  controller  performance  should  be  developed  based 
on  the  probability  of  various  atmospheric  perturbations  occurring.  This  method  may 
extend  further  to  include  the  probability  of  variations  in  entry  conditions  or  vehicle 
aerodynamic  characteristics. 

8)  A  new  performance  index  should  be  developed  which  will  minimize  AV  while  re¬ 
taining  a  level  of  robustness.  With  this  new  performance  index,  the  calculus  of  vari¬ 
ations  optimization  techniques  should  be  revisited  in  an  attempt  to  construct  a 
controller  which  computes  a  truly  optimal  solution. 
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APPENDIX  A 

IDEALIZED  MINIMUM  AV  OPTIMAL  SOLUTION 

A  numerical  gradient  technique  was  employed  to  determine  the  minimum  AV  solu¬ 
tion  for  a  nominal  Martian  aerobraking  maneuver'^.  The  MRSR  mission  scenario  calls 
for  the  aerobraking  maneuver  to  reduce  the  vehicle’s  velocity  relative  to  the  planet  using 
aerodynamic  drag  and  then  exit  the  atmosphere  on  an  elliptical  intermediate  orbit.  A  se¬ 
ries  of  propulsive  maneuvers  are  then  performed  to  transfer  the  vehicle  from  the  interme¬ 
diate  orbit  to  the  desired  final  orbit.  The  total  AV  required  to  transition  from  the 
intermediate  orbit  to  the  desired  orbit  is  determined  by  the  vehicle’s  atmospheric  exit  ve¬ 
locity  vector  and  is  a  good  measure  of  control  system  performance.  The  open  loop  solu¬ 
tion  presented  here  assumes  that  initial  conditions  as  well  as  all  pertinent  vehicle  and 
atmospheric  properties  are  known  precisely.  Limits  are  not  placed  on  trajectory  loads. 
Robustness  to  atmospheric  dispersions  is  not  considered  in  computing  this  optimal  solu¬ 
tion.  This  solution  produces  the  minimum  AV  attainable  to  transition  from  the  post  aero¬ 
braking  intermediate  orbit  to  the  desired  final  orbit  for  a  given  atmosphere,  vehicle  and 
entry  condition  and  is  used  as  a  benchmark  to  evaluate  the  performance  of  the  feedback 
controllers. 


Equations  of  Motion 

The  formulation  begins  with  the  equations  of  motion.  The  equations  of  motion  were 
presented  in  Chapter  111  but  are  repeated  again  here  for  completeness. 


dr  dh  . 

3?  =  J,  = 


(86) 
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dt 


-CoP-svJ 

2m 


(87) 


dt 


2mV 


cosO  - 


(88) 


Equation  (86)  is  simply  the  radial  velocity  in  terms  of  the  inertial  velocity  and  flight 
path  angle.  Equation  (87)  gives  the  time  rate  of  change  of  velocity  in  two  parts:  1)  the  ve¬ 
locity  loss  rate  due  to  aerodynamic  drag  and  2)  the  change  in  velocity  due  to  gravitational 
acceleration  (the  inertial  component).  Similarly,  Eq.  (88)  is  the  time  rate  of  change  in  the 
flight  path  angle  also  composed  of  two  parts:  1)  the  change  in  flight  path  angle  due  to  the 
component  of  aerodynamic  lift  in  the  vertical  plane  and  2)  the  change  in  flight  path  angle 
due  to  gravitational  acceleration  (the  inertial  component).  The  control  variable  O,  the 
bank  angle,  determines  the  amount  of  lift  exerted  in  the  vertical  plane  to  bend  the  trajecto¬ 
ry  and  change  the  flight  path  angle. 


Nondimensional  State  Variables 

Dimensionless  state  variables  are  introduced: 


/.//I, 

X  = 

^2 

= 

V/U\i/R) 

.^3 

y 

along  with  a  dimensionless  time  variable  x 

X  =  it/h^) 


(89) 


(90) 


The  equations  of  motion  may  now  be  written: 
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jfj  =  X2SinjC3 


(91) 


X2  =  -80x2  - 


{c-  1  +x,) 


^sinx^ 


(92) 


Aaxr 


^3  = 


-COS<I>  + 


COSXi 


C  -  1  +JC, 


^2  (c-l+x,)x 


(93) 


where  a  =  p/p^  =  exp  [  (-(/i-/iQ))/hS] ,  A  =  (p^Sh^Cj)  /  (2m) , 
B  =  (p^ShgC q)  /  (2m)  smA  c  =  R/h^. 


The  Performance  Index 

To  minimize  total  AV  required  to  transition  to  the  desired  orbit  it  is  sufficient  to  min¬ 
imize  the  exit  flight  path  angle  provided  the  apocenter  of  the  intermediate  orbit  equals  the 
desired  apocenter.  This  procedure  maximizes  the  pericenter  of  the  post-aero  orbit.  Two 
terminal  constraints  are  employed.  The  first  requires  the  final  altitude  to  be  the  atmospher¬ 
ic  interface  altitude  and  the  second  fixes  the  intermediate  orbit  apocenter.  The  cost  func¬ 
tion  is  therefore  the  exit  flight  path  angle  (J  =  yp  tmd  the  goal  is  to  minimize  the  cost 
function  subject  to 


V,  =  x^-  1  =  0 

and 


(94) 


/  a. 

V2  =  -  (^) 


2c 


■^2-  , 

^  c  -  1  +  Jt 


rc-  1  +JCi' 

R 

L  c 

X2(COSiX-^)^  =  0 


(95) 
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To  derive  Eq.  (95),  set  the  orbital  angular  momentum  at  exit  equal  to  the  angular  mo¬ 


mentum  at  apocenter. 


h  =  r^V^cosy^  =  r^V^. 


From  this  equation  solve  for  the  velocity  at  apocenter  in  terms  of  the  terminal  radius, 
velocity,  flight  path  angle  and  the  radius  of  apocenter. 


r^V^cosY_^ 


Equate  the  orbital  energy  at  exit  to  that  at  apocenter,  using  the  expression  for  veloci¬ 


ty  at  apocenter  from  above 


n  Vi 


2  r. 


Obtain  Eq.  (95)  after  some  algebra  and  after  replacing  the  physical  state  variables 


with  the  nondimensional  variables  given  in  Eqs.  (89)  and  (90). 


The  Numerical  Gradient  Technique 


This  problem  was  solved  using  a  first  order  numerical  gradient  procedure.  To  for¬ 
mulate  the  optimal  control  problem  begin  with  the  performance  index.  In  general  terms 


this  index  may  be  written 
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The  performance  index  is  then  augmented  with  penalty  functions  to  impose  the  ter¬ 
minal  constraints  and  the  equations  of  motion. 

Ja  =  I/*)  +  [Xl'T  [/•(/,«*)-/]}  Jr  (100) 

For  this  problem  [v]  is  a  2  x  1  column  matrix  of  constants,  {V(xy)  }  is  a  2  x  1 
column  matrix  of  terminal  constraints  given  by  Eqs.  (94)  and  (95)  and  [X]  is  a  3  x  1  time 
varying  matrix  of  Lagrange  multipliers  or  influence  functions,  {x}  =  {f{x,  u)  }  are  the 
3  first  order  differential  equations  of  motion  Eqs.  (91),  (92)  and  (93).  {«}  is  the  control 
variable  O.  To  customize  this  general  augmented  cost  function  for  the  problem  at  hand 
delete  L  (x*,  u*)  since  there  is  not  an  integral  term  in  our  performance  measure,  and  sub¬ 
stitute  7^  for  (J)  (xy)  to  obtain 

h  =  7/+  +||'{  -/))*  (101) 

The  numerical  solution  process  begins  with  a  guess  of  the  control  time  history.  The 
values  for  the  state  variables  are  computed  from  initial  conditions  and  then  integrated  for¬ 
ward  in  time  using  this  postulated  control  time  history.  Differential  equations  for  the 
Lagrange  multipliers,  which  are  developed  later,  are  used  to  integrate  [X]  backward  in 
time  beginning  with  the  value  of  [X]  computed  at  tj.  A  new  control  time  history  is  de¬ 
rived  by  setting  the  first  variation  in  the  augmented  cost  function  to  zero.  The  process  is 
repeated  until  the  terminal  constraints  are  satisfied  to  within  an  acceptable  tolerance. 


The  first  variation  in  is  formed 
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Integrate  [A,]  by  parts  to  obtain 

(103) 

Substitute  Eq.  (103)  into  Eq.  (102)  to  obtain 

'■'o "  ''''  %" 

Define  a  new  Lagrange  multiplier 

[X]T's[X^]'*'+[v]T[V]'^  (105) 

[X-^]  is  the  3x1  column  of  Lagrange  multipliers  normally  used  to  impose  the 
equations  of  motion  while  [A,']  is  dimensioned  3x2  and  contains  additions  to  the 
Lagrange  multipliers  which  arise  from  the  terminal  constraints.  The  first  variation  in 
may  now  be  written 


rt>y 


J)x 


^4  [VI  ^^-ix;i"-[v)Tix;, 


-r  T  cif  •  y  T  *  . ,  T 

+  [v]^[V]')^+[n  +[v]  [A.']  )5x 


+  (([Av-']'’'  +  rv]'*[X']'^)|^)6M]ryr 


(106) 


148 


The  performance  index  is  minimized  by  setting  the  first  variation  in  equal  to  zero. 
Choose  differential  equations  for  the  Lagrange  multipliers  so  that  the  coefficient  of  6jc 
goes  to  zero  to  obtain  the  two  equations 


and 


The  gradient  of  /  is 


where 


0 

sinxj . 

r2sinx 

a/2 

3/2 

9/2 

^/3 

5/3 

a/3 

dx^ 

al^ 

2csinx3 

^  —  + - r 

dxj  hS  (c_i+x,)3 


-2Bax2r 


0/2  CCOSAT^ 
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(107) 


(108) 


(109) 


(110) 


(111) 


(112) 
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8^2 


-2 


(c  1  +  Jf  j^)  X2^ 
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(120) 


Since  the  coefficients  of  5x  and  Sxj-  have  been  set  to  zero,  the  first  variation  of  re¬ 
duces  to 


=  87+  [v]'r{8\|r} 


(121) 


where 


57  = 


and 


(122) 


{5v}  =  J[(M^g^)5M]7r 


(123) 


Defining  two  new  variables 


(124) 


(125) 


is  a  scalar,  while  is  a  2  x  1  column  matrix  and 


^  =  jo  0  -Aax2sin<I^. 


(126) 
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To  make  as  small  as  possible,  choose  the  variation  in  the  control  5m  to  be 

bu  =  (127) 

K  \s  a  scalar  weight  which  fixes  the  relative  importance  placed  on  minimizing  the 
cost  function  versus  satisfying  the  terminal  constraints.  A  value  of  200  for  K  places  suffi¬ 
cient  weight  on  the  cost  function  and  still  allows  the  terminal  constraints  to  be  satisfied 
within  an  acceptable  tolerance.  By  substituting  Eqs.  (127)  and  (125)  into  Eq.  (123)  obtain 


{8v}  =  -Kjl' [A^l  [A^+  [V]’'A^]'^*.  (128) 

Again,  introduce  two  additional  variables 

{«}  =|;Ma^HA^1T*  (129) 

IQl-jlMA^llA^lT*  (130) 

Substitute  these  variables  into  Eq.  (128)  to  obtain 

{6v}  =  -K[g  +  Q[v]].  (131) 


To  drive  {8\}/}  to  zero,  choose  {5\jr}  =-{  ,  where  {yy}  is  the  value  of  the 

terminal  constraints,  computed  after  integrating  the  state  equations  forward,  solving  Eq. 
(131)  for  {V} 


(v)  = 


(132) 
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(135) 


Replace  <1)  with  and  \jr  with  the  expressions  given  in  Eqs.  (94)  and  (95)  to  obtain 


ay  av2  ^¥2 


(136) 


Use  the  new  control  time  history  {Eq.  (133)),  along  with  the  change  in  {Eq. 
(136)},  to  again  integrate  the  stale  equations  forward.  Compute  the  terminal  value  of  the 
Lagrange  multipliers  and  integrate  these  backwards  in  time,  then  recompute  the  control 
time  his.  iry.  This  process  is  repeated  until  the  terminal  constraints  are  satisfied  within  an 
acceptable  error  bound.  The  final  apocenter  altitude  was  required  to  be  within  5  nm  of  the 
target  va;ue  while  the  terminal  altitude  was  required  to  be  within  25,000  ft  of  the  defined 
atmospheric  interface  altitude.  Apocenter  errors  of  5  nm  require  very  little  AV  to  correct 
and  are  attainable  using  this  optimization  method  although  thirty  or  more  iterations  may 
be  required  to  converge  this  closely.  The  25,000  ft  terminal  altitude  error  band  was  chosen 
because  the  aerodynamic  effects  decrease  exponentially  with  altitude  and  are  almost  negli- 
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gible  at  altitudes  above  300,000  ft.  To  converge  closer  than  25,000  ft  to  the  selected  atmo¬ 
spheric  interface  altitude  of  410,105  ft  (125  km)  requires  many  more  iterations. 
Furthermore,  the  terminal  altitude  generally  converges  toward  the  target  altitude  from 
above. 


Conjugate  Gradient  Projection  Method 

The  conjugate  gradient  projection  method  was  employed  to  speed  convergence  of 
this  problem'*®’  '**.  The  gradient  obtained  in  Eq.  (127)  was  again  used  in  this  method  to 
compute  the  search  direction  for  correcting  the  control  variable.  However,  after  the  first 
control  update  the  previous  search  direction  is  used  in  conjunction  with  the  computed  gra¬ 
dient  to  give  the  problem  near  second-order  convergence  characteristics.  The  procedure 
follows.  First,  compute  the  gradient  direction  using  Eq.  (127) 

Si  =  -5m  (137) 

where  the  i  subscript  denotes  the  /th  iteration  of  control  updates.  Next,  compute  the  search 
direction 


=  -Si 


(*38) 


for  the  first  iteration,  while  for  subsequent  iterations 


{gj,  gi) 


(139) 


where  {a,  b)  is  the  inner  product  of  a  and  b. 

Once  the  search  direction  is  determined,  it  is  necessary  to  properly  scale  the  magni¬ 


tude  of  the  correction. 
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t“new>  =  ^“old> 

A  line  search  employing  a  Newton  iterative  scheme  is  used  to  scale  k  so  that  ^2 
minimized.  A  value  for  A:  is  chosen  (Jtj)  and  the  equations  of  motion  are  integrated  for¬ 
ward  and  the  values  of  the  terminal  constraints  are  determined.  Then,  a  larger  value  of  k  is 
chosen  {k^  and  the  equations  of  motion  are  again  integrated  forward  and  the  terminal 
constraint  is  determined.  A  new  value  for  k  is  selected  using 


A:  .  ,  =  kj - - - ^ - 

7  +  1  7 


V. 


(141) 


This  iteration  is  repeated  until  the  terminal  constraint  \\f2  is  within  an  acceptable  tol¬ 
erance  of  zero.  This  tolerance  is  again  computed  by  requiring  the  terminal  apocenter  error 
to  be  within  5  nm  of  the  target.  A  new  gradient  is  then  computed  and  a  new  search  direc¬ 
tion  is  found.  The  line  search  is  then  repeated. 


The  conjugate  gradient  procedure  requires  computation  of  the  gradient  as  before, 
then  computation  of  the  search  direction  using  Eqs.  (138)  or  (139).  Finally,  a  line  search  is 
employed  to  determine  the  desired  magnitude  of  the  correction.  This  procedure  is  repeat¬ 
ed  until  the  inner  product  of  the  computed  gradient  is  sufficiently  small  (less  than  10'^). 


Both  the  conjugate  gradient  technique  and  the  numerical  gradient  technique  pro¬ 
duced  acceptable  results  for  determining  optimum  (minimum  AV)  performance,  though 
the  conjugate  gradient  procedure  converged  somewhat  faster  than  the  numerical  gradient 
technique.  The  conjugate  gradient  method  required  computation  of  the  gradient  generally 
only  four  or  five  times.  The  line  search,  however,  was  a  slow  expensive  process  requiring 
as  many  as  ten  iterations  for  each  search  direction.  Overall,  the  conjugate  gradient  tech¬ 
nique  did  converge  faster  than  the  numerical  gradient  technique  but  only  by  about  50%. 
To  be  used  to  compute  optimal  control  time  histories  onboard  the  satellite  in  real  time,  a 
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solution  must  be  obtained  at  least  two  orders  of  magnitude  faster  than  either  the  numerical 
gradient  or  conjugate  gradient  techniques  currently  achieve. 
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APPENDIX  B 

ANALYTIC  PREDICTOR  CORRECTOR  DERIVATIONS 

An  Analytic  Predictor/Corrector  (APC)  algorithm  for  the  MRSR  was  adapted  from 
the  Aeroassist  Flight  Experiment  (AFE)  controller^*  This  controller  has  two  phases. 
There  is  an  equilibrium  glide  phase  during  the  first  part  of  the  trajectory,  where  aerody¬ 
namic  loads  are  the  primary  concern.  At  a  predetermined  velocity  the  controller  switches 
to  a  predictor/corrector  algorithm  for  the  exit  phase.  The  predictor  step  assumes  constant 
altitude  rate  and  analytically  integrates  the  trajectory  forward.  Then,  it  corrects  for  the  fi¬ 
nal  phases  of  the  trajectory  where  constant  altitude  rate  cannot  be  maintained.  The  correc¬ 
tor  step  adjusts  altitude  rate  to  target  the  desired  apoapsis,  thereby  minimizing  AV  required 
for  insertion  into  the  desired  low  Mars  orbit.  Also  included  in  this  Appendix  is  the  devel¬ 
opment  of  Fitzgerald’s  Hybrid  Predictor  Corrector*  *  exit  phase  which  uses  density  values 
derived  from  accelerometer  measurements  during  entry  to  obtain  a  better  estimate  of  the 
velocity  loss  due  to  aerodynamic  drag  during  the  exit  phase. 


Equilibrium  Glide  Phase 


The  equilibrium  glide  phase  of  the  APC  controller  seeks  an  equilibrium  condition 
with  the  vehicle  following  a  reference  dynamic  pressure  path.  Equilibrium  is  established 
by  requiring  the  in-plane  portion  of  lift  to  balance  inertial  and  gravitational  forces. 


C^^cos<I)  =  wcosy- 


(142) 


For  small  flight  path  angles,  the  radius  of  curvature  approximates  the  vehicle’s 
orbital  radius.  Using  this  approximation,  the  bank  angle  required  to  maintain  equilibrium 
is 


cosO  = 


(143) 


The  reference  dynamic  pressure  is  calculated  as  a  multiple  K-  of  the  dynamic  pres¬ 
sure  required  to  maintain  equilibrium  with  the  lift  vector  oriented  down.  Previous 
works^’^^  have  recommended  Rg  =  1-33  for  an  Earth  aerobraking  vehicle.  However,  a 
value  of  4.5  for  K-  provides  additional  robustness  for  the  MPC.  Additional  discussion  of 

fi 

this  choice  may  be  found  in  Chapter  II  on  page  32. 


^ref 


(144) 


To  prevent  overshoot,  an  altitude  damper  is  included  in  the  control  equation  giving 
the  following  commanded  bank  angle: 


<D 
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acos 
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-G.h  +  G-{q-q^f) 


(145) 


Exit  Phase 

After  the  vehicle  slows  to  a  predetermined  velocity,  the  controller  transitions  to  a 


predictor/corrector  algorithm  which  targets  the  desired  apocenter  following  atmospheric 
exit.  The  analytical  relationship  used  to  compute  velocity  loss  during  the  remainder  of  the 
aerodynamic  phase  is  developed  by  considering  only  the  aerodynamic  drag 
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dV^  O.SpVjSC^ 

~dr  ""  m  * 

0.55C^ 

Lumping  constants  C  =  - and  assuming  an  exponential  atmosphere. 


^  =  -cv^D 

dt 


Rearranging  terms  and  replacing  dt  with  ^  we  get 

h 


dV^ 


-(/i-V/hS^ 
h  ' 


(146) 
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The  open  loop  optimal  minimum  AV  solutions  computed  using  the  methods  of  Ap¬ 
pendix  A  show  that  after  passing  the  pericenter,  h  is  very  nearly  a  constant.  By  assuming 
h  is  constant,  may  be  determined  at  any  future  altitude. 


-(h^-hQ)/hS  -(/i-/»o)/hS^ 

'(e  "  -e  ) 


-.-1 


(149) 


A  correction  must  be  added  to  account  for  the  kinetic-potential  energy  interchange. 

=  i‘'?  +  2(4'A-8A)-''/  050) 

Since  orbital  calculations  rely  on  inertial  velocity  rather  than  relative  velocity,  it  is 
desirable  to  compute  the  exit  inertial  velocity.  A  good  approximation  is  that  the  inertial 
and  relative  velocity  differ  by  a  constant  throughout  the  trajectory  and  the  correction  is 
simply 


(151) 
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The  inertial  exit  velocity  can  be  expressed  as 


^rx  (k/p)  '^^^x il/r)  * 


(152) 


As  altitude  increases  and  density  decreases  aerodynamics  become  less  dominant  and 
the  trajectory  becomes  more  an  inertial  orbital  trajectory  following  a  Keplerian  elliptical 
orbit.  Soon,  it  becomes  impossible  to  maintain  a  constant  altitude  rate;  h  will  increase 
even  with  full  lift  down.  The  APC  algorithm  assumes  that  at  some  predetermined  switch 
altitude  h  begins  increasing  linearly  until  atmospheric  exit.  The  altitude  acceleration  at 
exit  can  be  computed  by  summing  forces  in  the  vertical  direction  at  exit  assuming  the  lift 
vector  is  oriented  down. 


K  =  'Sx-  (0.5p/?/C^) /m  (153) 

Assuming  h  increases  linearly  with  time  after  the  switch  altitude,  a  quadratic  in  At  is 
written  for  the  altitude  during  the  final  segment  of  the  trajectory  as  the  vehicle  nears  atmo¬ 
spheric  exit. 

+  +  (154) 

Solving  for  At  and  multiplying  by  h  we  get  the  change  in  altitude  rate. 

Alt  ~  J^const  ^'^^switch  ^const  (155) 

The  altitude  rate  at  exit  is  the  constant  altitude  rate  plus  the  change  in  altitude  rale. 

/ij.  —  ftfonst  (156) 


With  this  altitude  rate  and  the  exit  velocity  calculated  in  Eq.  (152)  the  apocenter  of 
the  orbit  following  exit  may  be  calculated.  By  iterating  h  and  predicting  the  apocenter,  an 
aporopriate  h  may  be  found  which  should  yield  the  desired  post-aerobraking  orbit  apo¬ 


center. 


To  gain  a  measure  of  robustness  for  both  density  deviations  and  deviations  in  Cp,  a 
density  filter  is  incorporated  which  uses  density  derived  from  measured  drag  deceleration 


Pd  = 


2Vm 

CpSVj 


(157) 


This  derived  density  is  divided  by  the  density  predicted  for  the  current  altitude  using 
a  standard  exponential  atmosphere.  The  result  is  filtered  using  a  low  pass  filter  to  remove 
high  frequency  density  deviations  which  would  have  minimal  effect  on  the  post-aerobrak¬ 
ing  apocenter. 

^p=  (*-^)^p  +  ^(P/Pn.odel)  (158) 

The  resulting  filtered  density  multiplier  is  multiplied  by  pQ  during  the  predictor  step. 
As  noted  by  Gamble,  et  aP‘^,  will  compensate  for  unce.dainties  in  the  aerodynamic 
drag  coefficient  Cp  as  well  as  density  uncertainties.” 


The  control  equation  for  the  exit  phase  is  very  similar  to  the  equilibrium  glide  con¬ 
trol  equation.  The  major  differences  are  the  inclusion  of  a  desired  altitude  rate  instead  of 
simply  an  altitude  damper  and  the  elimination  of  the  reference  dynamic  pressure  term. 
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Hybrid  Predictor  Corrector 

Hybrid  Predictor  Corrector  is  a  name  first  coined  by  Fitzgerald**  to  describe  an  im¬ 
provement  introduced  in  the  exit  phase  of  the  original  APC  algorithm.  Fitzgerald’s  ap¬ 
proach  is  to  derive  density  from  accelerometer  measurements  at  discrete  intervals  during 
the  descent  into  the  atmosphere.  The  exit  phase  then  fits  an  exponential  density  curve  to 
each  altitude  band  between  the  discrete  points.  The  velocity  loss  through  the  atmosphere 
due  to  aerodynamic  drag  is  then  calculated  as  the  summation  of  the  velocity  loss  through 
each  altitude  band. 


The  rationale  for  this  improvement  is  that  accelerometer-generated  density  measure¬ 
ments  taken  during  the  entry  phase  of  the  aerobraking  maneuver  are,  quite  likely,  a  reason¬ 
able  estimate  of  the  atmospheric  density  function  available  for  the  exit  phase  of  the 
trajectory.  These  measurements  will  be  close,  in  both  space  and  time,  to  the  exit  phase  of 
the  flight  and  will  hopefully  produce  a  good  estimate  of  the  density  to  be  encountered  dur¬ 
ing  the  exit  phase.  In  this  development  p  |  is  the  density  which  was  measured  at  the  lower 
edge  of  the  current  altitude  band,  is  the  altitude  at  which  this  measurement  was  taken. 
P2  is  the  density  which  was  measured  at  the  upper  edge  of  this  altitude  band  at  altitude  /j2. 
hS  is  the  scale  height  for  the  atmosphere  band  computed  between  the  two  density  mea¬ 
surements. 


hS  = 


logCp^/p,) 

//,  -/j2 


-1 


(160) 


To  use  this  modified  atmosphere  in  the  predictor  step,  rewrite  Eq.  (148) 


dV. 


,2 


-(/i-*o)/hS^ 

ii 


(161) 


This  equation  may  be  integrated  assuming  a  constant  altitude  rate  to  give  the  velocity 
loss  due  to  atmospheric  drag  between  two  arbitrary  altitudes  /j  j  and  hj- 
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''r2  = 


1  ("<^PohS^^^-(A2-/,jj)/hs  _^-(y.,-/.o)/hs^ 


-1 


(162) 


The  velocity  at  altitude  /12  in  terms  of  the  velocity  at  and  the  density  at  each  loca¬ 
tion  is  defined  by 


^.2  = 


L^rl 


1  f^PlhS^^^_(|,2-A,)/hS  _^-(/i,-li,)/hS^ 


1-1 
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and,  with  the  scale  height  as  previously  calculated 


^^^2  = 


1  Cp,(Aj-/l2) 


y^r\  Alogfp^/pj)  iPl 


1-1 


(164) 


This  equation  gives  the  relative  velocity  at  ^2  ns  a  function  of  the  relative  velocity  at 
/ij  and  the  densities  and  altitudes  at  the  two  locations.  The  method  for  employing  this  fea¬ 
ture  in  the  predictor  step  of  the  control  algorithm  is  to  first  use  the  velocity,  density  and  al¬ 
titude  at  the  current  satellite  location  as  the  subscript  1  variables  and  to  predict  the  velocity 
at  the  next  interval  where  density  measurements  were  stored  during  the  entry  using  that  al¬ 
titude  and  that  density  as  the  subscript  2  variables.  Then,  that  velocity  may  be  used  to 
compute  the  velocity  at  the  next  altitude  band  using  the  lower  stored  density  and  altitude 
values  as  subscript  1  variables  and  the  next  higher  density  and  altitude  measurements  as 
subscript  two  variables.  This  procedure  is  repeated  until  the  exit  relative  velocity  is  com¬ 
puted;  that  exit  velocity  is  then  handled  exactly  as  it  was  for  the  APC. 


APPENDIX  C 


ENERGY  CONTROLLER  DERIVATION 


The  Energy  Controller^'*  was  also  adapted  to  control  the  MRSR  vehicle.  The  Energy 
Controller  defines  a  new  variable,  the  energy  gain,  as  the  ratio  of  energy  rate  to  energy  er¬ 
ror.  The  energy  gain  is  controlled  so  that  Keplerian  energy  approaches  the  commanded 
value  as  energy  rate  goes  to  zero  at  atmospheric  exit,  directing  the  vehicle  orbit  to  the  de¬ 
sired  apocenter.  An  analytic  relationship  is  used  to  convert  energy  gain  into  an  altitude 
rate  command.  The  altitude  rate  error  is  used  to  compute  a  desired  altitude  acceleration 
which  leads  to  a  desired  bank  angle. 

Development  of  this  controller  begins  with  the  energy  to  mass  ratio. 


2  r 

tl  ^ 


(165) 


To  calculate  the  commanded  energy  at  exit  the  conical  equivalent  is  substituted  for 


Rp  =  a{\'e) 


(166) 


so  that  the  desired  energy  at  exit  may  be  computed  from 
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Energy  Gain 


The  main  control  variable,  energy  gain  E  ,  is  defined  as  the  ratio  of  E  to  energy  er- 

o 


ror. 


= 


(168) 


The  energy  gain  is  controlled  so  that  E  approaches  the  commanded  value  exponen¬ 
tially  as  the  energy  rate  approaches  zero. 


~  ^gcO'^Egcxf 


(169) 


EgcO  initial  energy  gain  command  and  Eg^.^  is  the  desired  energy  gain  rate  at 
atmospheric  exit.  Eg^.^  will  be  derived  later.  During  the  atmospheric  entry  phase  a  first- 
order  controller  is  u.sed  for  the  commanded  energy  gain  rate. 


^gc  "  Eg^^  +  k^(E^^  Eg) 


(170) 


Later  the  proportional  term  is  dropped  so  that 


Egc  Eg^^. 


(171) 


DV  I  ^ 

Differentiating  Eq.  (168),  E  =  — V'  =  J2E — D  =  - ^ 

.  A/i 

p  =  p^e  results  in  the  following: 
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V  _  E  \ih 

y  yp-  (;.v/)2 

(174) 

D  p 

(175) 

1 

II 

•Q.IQ. 

(176) 

which  can  be  solved  for  h. 

.  /jS[3(£/v2)+£:  -V^pl 

h  —  s  s  g 

(rV)2 

(177) 

3nhS 

Since - :r  is  small  compared  to  1,  drop  this  term  and  substitute  £„r  for  £„ 

(rVO^  ^  ^ 

tain  the  altitude  rate  command. 

to  ob- 

(178) 

Altitude  Acceleration  Command 

The  commanded  altitude  acceleration  is  a  first-order  control  on  the  altitude  rate  error 

plus  a  lead  term. 

h^.  =  hi  +  0.0S(h^-  h) 

(179) 
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The  lead  term  is  obtained  by  differentiating  Eq.  (178)  obtaining 


dh 


c  .  ^  E  .  ^  .  ^gc^g  ^^gc 


_ 

2  dt  i 


(180) 


Using  Eq.  (172),  assuming  Eg  -  Eg^,  differentiating  Eq.  (170)  with  Eg^^  =  0,  and 
assuming  Eg^.  =  Eg^j^  gives  the  lead  term  for  the  commanded  acceleration  equation. 


hi  =  hS 


(181) 


Differentiating  fi  gives  E 
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h  =  hS 


UoJ  ~  e:, 

^  g. 
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At  exit  E  and  Eg  approach  zero;  so  Eq.  (182)  gives  the  exit  altitude  acceleration, 
which  also  equates  to  the  orbital  dynamics  expression  for  altitude  acceleration  at  exit 


hex  =  hS 


-gcx 


+  E 


gcx 


V^^cos^Y  ^ 
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By  solving  eq.  (165)  for  V 


ex 


Vt.  =  2E  + 


2\i 


ex 


(184) 


and  substituting  this  result  into  the  previous  equation  to  obtain 


hex  ~  +  • 


(185) 
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Now,  equating  from  above  and  Eq.  (183)  and  solving  for 


E  = 

^gcx 


|()iS)2  +  4/j5^ 

(186) 


The  total  vertical  acceleration  may  be  written  as 


(Vcosy)  |i  LcosOcosy  Dsiny 

h  - - ^  H - 

r  m  m 
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Bank  Angle  Command 

Solve  for  the  commanded  bank  angle  by  substituting  the  commanded  altitude  accel¬ 
eration  in  Eq.  (187). 


cos<I>  = 
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(188) 
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